R Workflow

Author

Frank Harrell

Published

May 18, 2022

This article outlines analysis project workflow that I’ve found to be efficient in making reproducible research reports using R with Rmarkdown and now Quarto. I start by covering the creation of annotated analysis files and running descriptive statistics on them with goals of understanding the data and the quality and completeness of the data. Functions in the Hmisc package are used to annotate data frames and data tables with labels and units of measurement and to produce tabular and graphical statistical summaries. Several examples of processing and manipulating data using the data.table package are given. Finally, examples of caching results and doing parallel processing are presented.

Comments

Assignment Operator

You assign an R object to a value using the assignment operator <- or the equal sign. <- is read as “gets”.

x <- y
d <- read.csv('mydata.csv')
x = y

Object Types

Everything in R is an object. Some primitive types of objects in R are below.

Type Meaning
integer whole numbers
logical values of TRUE or FALSE
double floating point non-whole numbers
character character strings
function code defining a function

In the table below, objects of different shapes are described. rows and cols refers to vectors of integers or logicals, or if the elements of the object are named, character strings.

Type Example Values Retrieved By
scalar x <- 3 x
vector y <- c(1, 2, 5) y[2] (2), y[2:3] (2, 5), y[-1] (2, 5), y[c(TRUE,FALSE,TRUE)] (1)
named vector y <- c(a=1, b=2, d=5) y[2] (2), y['b'] (2), y[c('a','b')] (1, 2)
matrix y <- cbind(1:3, 4:5) y[rows,cols], y[rows,] (all cols), y[,cols] (all rows)
list x <- list(a='cat', b=c(1,3,7)) x$a (‘cat’), x[[1]] (‘cat’), x[['a']] (‘cat’)

Named vectors provide an extremely quick table lookup and recoding capability.

list objects are arbitrary trees and can have elements nested to any level. You can have lists of lists or lists of data frames/tables.

Vectors can be of many different types when a class is added to them. Two of the most common are Dates and factors. Character strings are handled very efficiently in R so there is not always a need to store categorical variables as factors. But there is one reason: to order levels, i.e., distinct variable values, so that tabular and graphical output will list values in a more logical order than alphabetic. A factor variable has a levels attribute added to it to accomplish this. An example is x <- factor(x, 1:3, c('cat', 'dog', 'fox')) where the second argument 1:3 is the vector of possible numeric values x currently takes on (in order) and the three character strings are the corresponding levels. Internally factors are coded as integers, but they print as character strings.

Rectangular data objects, i.e., when the number of rows is the same for every column (variable), can be represented by matrices, data.frames, and data.tables. In a matrix, every value is of the same type. A data.frame or a data.table is an R list that can have mixtures of numeric, character, factor, dates, and other object types. A data.table is also a data.frame but the converse isn’t true. data.tables are handled by the R data.table package and don’t have row names but can be indexed, are much faster to process, and have a host of methods implemented for aggregation and other operations. data.frames are handled by base R.

Subscripting

Examples of subscripting are given above. Subscripting via placement of [] after an object name is used for subsetting, and occasionally for using some elements more than once:

x <- c('cat', 'dog', 'fox')
x[2:3]
[1] "dog" "fox"
x[c(1, 1, 3, 3, 2)]
[1] "cat" "cat" "fox" "fox" "dog"

Subscripting a variable or a data frame/table by a vector of TRUE/FALSE values is a very powerful feature of R. This is used to obtain elements satisfying one or more conditions:

x <- c(1, 2, 3, 2, 1, 4, 7)
y <- c(1, 8, 2, 3, 8, 9, 2)
x[y > 7]
[1] 2 1 4

The last line of code can be read as “values of x such that y > 7”.

Branching and If/Then

Decisions Base on One Scalar Value

Common approaches to this problem are if and switch.

type <- 'semiparametric'
f <- switch(parametric     = ols(y ~ x),
            semiparametric = orm(y ~ x),
            nonparametric  = rcorr(x, y, type='spearman'),
            { z <- y / x
              c(median=median(z), gmean=exp(mean(log(z))) } )
# The last 2 lines are executed for any type other than the 3 listed
f <- if(type == 'parametric')    ols(y ~ x)
  else
    if(type == 'semiparametric') orm(y ~ x)
  else
    if(type == 'nonparametric')  rcorr(x, y, type='spearman')
  else {
    z <- y / z
    c(median=median(z), gmean=exp(mean(log(z)))
  }

What is inside if( ) must be a single scalar element that is evaluated to whether it’s TRUE or FALSE.

Series of Separate Decisions Over a Vector of Values

The ifelse or data.table::fifelse functions are most often used for this, but data.table::fcase is a little better. Here’s an example.

x <- c('cat', 'dog', 'giraffe', 'elephant')
type <- ifelse(x %in% c('cat', 'dog'), 'domestic', 'wild')
type
[1] "domestic" "domestic" "wild"     "wild"    
require(data.table)
fcase(x %in% c('cat', 'dog'), 'domestic', default='wild')
[1] "domestic" "domestic" "wild"     "wild"

Functions

Even new R users can benefit from writing functions to reduce repetitive coding. A function has arguments and these can have default values for when the argument is not specified by the user when the function is called. Here are some examples. One line functions do not need to have their bodies enclosed in {}.

cuberoot <- function(x) x ^ (1/3)
cuberoot(8)
[1] 2
g <- function(x, power=2) {
  u <- abs(x - 0.5)
  u / (1. + u ^ power)
}
g(3, power=2)
[1] 0.3448276
g(3)
[1] 0.3448276
# Function to make mean() drop missing values without our telling it
mn <- function(x) mean(x, na.rm=TRUE)
# Function to be used throughout the report to round fractional values 
# by a default amount (here round to 0.001)
rnd <- function(x) round(x, 3)
# edit the 3 the change rounding anywhere in the report
# The following simple function saves coding when you need to recode multiple
# variables from 0/1 to no/yes.
yn <- function(x) factor(x, 0:1, c('no', 'yes'))

Report Formatting

A state-of-the-art way to make reproducible reports is to use a statistical computing language such as R and its knitr package in conjunction with either RMarkdown or Quarto, with the latter likely to replace the former. Both of the report-making systems allow one to produce reports in a variety of formats including html, pdf, and Word. Html is recommended because pages can be automatically resized to allow optimum viewing on devices of most sizes, and because html allows for interactive graphics and other interactive components. Pdf is produced by converting RMarkdown or Quarto-produced markdown elements to \(\LaTeX\).

Report formatting is very much enhanced by using variable attributes such as labels and units of measurement that are not considered in base R. Methods for better annotating output using labels and units are given below.

This document can serve as a template for using R with Quarto; one can see the raw script by clicking on Code at the top right of the report. When one has only one output format target, things are fairly straightforward except for the following. To make use of specialized functions that produce html or \(\LaTeX\) markup, one often has to put results='asis' in the code chunk header to keep the system from disturbing the generated html or \(\LaTeX\) markup so that it will be typeset correctly in the final document. This process works smoothly but creates one complication: if you print an object in the same code chunk that produces plain text, the system will try to typeset it in html or \(\LaTeX\). To prevent this from happening you either need to split the chunks into multiple chunks (some with results='asis' and some not) or you need to make it clear that parts of the output are to be typeset verbatim. To do that a simple function pr can sense if results='asis' is in effect for the current chunk. If so, the object is surrounded by the markdown verbatim indicator—three consecutive back ticks. If not the object is left alone. pr is defined in the marksupSpecs$markdown$pr object, so you can bring it to your session by copying into a local function pr as shown below, which has a chunk option results='asis' to show that verbatim output appears anyway.

If the argument obj to pr is a data frame or data table, variables will be rounded to the value given in the argument dec (default dec=3) before printing. If you specify inline=x the object x is printed with cat() instead of print(). inline is more for printing character strings.

An example of something that will not render correctly due to results='asis' being in the chunk header (needed for html(...)):

html(describe(d))
table(sex, country)
# use pr(obj=table(sex, country))
# to fix
require(Hmisc)
pr <- markupSpecs$markdown$pr
x <- (1:5)/7
pr('x:', x)

x: 

[1] 0.1428571 0.2857143 0.4285714 0.5714286 0.7142857
pr(obj=x)
[1] 0.1428571 0.2857143 0.4285714 0.5714286 0.7142857
pr(inline=paste(round(x,3), collapse=', '))

0.143, 0.286, 0.429, 0.571, 0.714

Instead of working to keep certain outputs verbatim you can use knitr::kable() to convert verbatim output to markdown.

Multi-Output Format Reports

To allow one report to be used to render multiple output formats, especially html and pdf, it is helpful to be able to sense which output format is currently in play, and to use different functions or options to render output explicitly for the current format. Here is how to create variables that can be referenced simply in code throughout the report, and to invoke the plotly graphics package if output is in html to allow interactivity. A small function ggp is defined so that if you run any ggplot2 output through it, the result will be automatically converted to plotly using the ggplotly function, otherwise it is left at standard static ggplot2 output if html is not the output target.

See this for examples of articles rendered in both html and PDF
outfmt <- if(knitr::is_html_output ()) 'html'  else 'pdf'
markup <- if(knitr::is_latex_output()) 'latex' else 'html'
ishtml <- outfmt == 'html'
if(ishtml) require(plotly)
ggp <- if(ishtml) ggplotlyr else function(ggobject, ...) ggobject
# See below for more about ggplotlyr (a front end for ggplotly that can
# correct a formatting issue with hover text)

The Hmisc, rms, and rmsb packages have a good deal of support for creating \(\LaTeX\) output in addition to html. They require some special \(\LaTeX\) packages to be accessed. In addition, if using any of Quarto’s nice features for making marginal notes, there is another \(\LaTeX\) package to attach. Below you’ll find what needs to be added to the yaml prologue at the top of your script if using Quarto. You have to modify pdf-engine to suit your needs. I use luatex because it handles special unicode characters. In the future (approximately July 2022) a bug in Pandoc will be fixed and you can put links-as-notes: true in the yaml header instead of redefining href and linking in hyperref.

format:
  html:
    self-contained: true
    . . .
  pdf:
    pdf-engine: lualatex
    toc: false
    number-sections: true
    number-depth: 2
    top-level-division: section
    reference-location: document
    listings: false
    header-includes:
      \usepackage{marginnote, here, relsize, needspace, setspace, hyperref}
      \renewcommand{\href}[2]{#2\footnote{\url{#1}}}

The href redefinition above turns URLs into footnotes if running \(\LaTeX\).

There is one output element provided by Quarto that will not render correctly to \(\LaTeX\): a marginal note using the markup .column-margin. To automatically use an alternate in-body format, define a function that can be used for both typesetting formats.

mNote <- if(ishtml) '.column-margin'
  else
                    '.callout-note appearance="minimal"'

Then use r mNote enclosed in back ticks in place of the .column-margin callout for generality.

As done with various Hmisc and rms package functions, one can capitalize on Hmisc’s special formatting of variable labels and units when constructing tables in \(\LaTeX\) or html. The basic constructs are shown in the code below.

# Retrieve a set of markup functions depending on typesetting format
# See below for definition of ishtml
specs    <- markupSpecs[[if(ishtml) 'html' else 'latex']]
# Hmisc markupSpecs functions create plain text, html, latex,
# markdown, or plotmath code
varlabel <- specs$varlabel  # retrieve an individual function
# Format text describing variable named x
# hfill=TRUE typesets units to be right-justified in label
# Use the following character string as a row label
# Default specifies the string to use if there is not label
# (usually taken as the variable name)
varlabel(label(x, default='x'), units(x), hfill=TRUE)

File Directory Structure

I have a directory for each major project, and put everything in that directory (including data files) except for graphics files for figures, which are placed in their own subdirectory underneath the project folder. The directory name for the project includes key identifying information, and files within that directory do not contain project information in their names, nor do they contain dates, unless I want to freeze an old version of an analysis script.

With Quarto I specify that html files are to be self-contained, so there are no separate graphics files.

For multi-analyst projects or ones in which you want to capture the entire code history, having the project on github is worthwhile.

Analysis File Creation

I typically create a compact analysis file in a separate R script called create.r and have it produce compressed R binary data.frame or data.table .rds files using saveRDS(name, 'name.rds', compress='xz'). Then I have an analysis script named for example a.qmd (for Quarto reports) or a.Rmd (for RMarkdown reports) that starts with d <- readRDS('name.rds').

Templates for analysis reports are here, and a comprehensive report example may be found here.

When variables need to be recoded, have labels added or changed, or have units of measurement added, I specify those using the Hmisc package upData function.

Variable labels and units of measurement are used in special ways in my R packages. This will show up in the describe and contents function outputs below and in axis labels for graphics.

To facilitate some operations requiring variable names to be quoted, define a function .q to quote them automatically. .q is like the Hmisc function Cs but also allows elements to be named. It will be in Hmisc 4.7-1.

.q <- function(...) {
  s <- sys.call()[-1]
  w <- as.character(s)
  n <- names(s)
  if(length(n)) names(w) <- n
  w
}

.q(a, b, c, 'this and that')
[1] "a"             "b"             "c"             "this and that"
.q(dog=a, giraffe=b, cat=c)
    dog giraffe     cat 
    "a"     "b"     "c"

Here is an upData example:

# Function to recode from atypical coding for yes/no in raw data
yn <- function(x) factor(x, 0:1, c('yes', 'no'))
d <-
  upData(d,
         rename = .q(gender=sex, any.event=anyEvent),
         posSE    = yn(posSE),
         newMI    = yn(newMI),
         newPTCA  = yn(newPTCA),
         newCABG  = yn(newCABG),
         death    = yn(death),
         hxofHT   = yn(hxofHT),
         hxofDM   = yn(hxofDM),
         hxofCig  = factor(hxofCig, c(0, 0.5, 1),
                           c('heavy', 'moderate', 'non-smoker')), 
         hxofMI   = yn(hxofMI),
         hxofPTCA = yn(hxofPTCA),
         hxofCABG = yn(hxofCABG),
         chestpain= yn(chestpain),
         anyEvent = yn(anyEvent),
         drop=.q(event.no, phat, mics, deltaEF,
                 newpkmphr, gdpkmphr, gdmaxmphr, gddpeakdp, gdmaxdp,
                 hardness),
         labels=c(
           bhr       = 'Basal heart rate',
           basebp    = 'Basal blood pressure',
           basedp    = 'Basal Double Product bhr*basebp',
           age       = 'Age',
           pkhr      = 'Peak heart rate',
           sbp       = 'Systolic blood pressure',
           dp        = 'Double product pkhr*sbp',
           dose      = 'Dose of dobutamine given',
           maxhr     = 'Maximum heart rate',
           pctMphr   = 'Percent maximum predicted heart rate achieved',
           mbp       = 'Maximum blood pressure',
           dpmaxdo   = 'Double product on max dobutamine dose',
           dobdose   = 'Dobutamine dose at max double product',
           baseEF    = 'Baseline cardiac ejection fraction',
           dobEF     = 'Ejection fraction on dobutamine', 
           chestpain = 'Chest pain', 
           ecg       = 'Baseline electrocardiogram diagnosis',
           restwma   = 'Resting wall motion abnormality on echocardiogram', 
           posSE     = 'Positive stress echocardiogram',
           newMI     = 'New myocardial infarction',
           newPTCA   = 'Recent angioplasty',
           newCABG   = 'Recent bypass surgery', 
           hxofHT    = 'History of hypertension', 
           hxofDM    = 'History of diabetes',
           hxofMI    = 'History of myocardial infarction',
           hxofCig   = 'History of smoking',
           hxofPTCA  = 'History of angioplasty',
           hxofCABG  = 'History of coronary artery bypass surgery',
           anyEvent  = 'Death, newMI, newPTCA, or newCABG'),
         units=.q(age=years, bhr=bpm, basebp=mmHg, basedp='bpm*mmHg',
           pkhr=mmHg, sbp=mmHg, dp='bpm*mmHg', maxhr=bpm,
           mbp=mmHg, dpmaxdo='bpm*mmHg', baseEF='%', dobEF='%',
           pctMphr='%', dose=mg, dobdose=mg)
         )

saveRDS(d, 'stressEcho.rds', compress='xz')

# Note that we could have automatically recoded all 0:1 variables
# if they were all to be treated identically:

for(x in names(d)) 
  if(all(d[[x]] %in% c(0,1))) d[[x]] <- yn(d[[x]])

Sometimes metadata comes from a separate source. Suppose you imported a data frame d but have also imported a data frame m containing metadata: the same variable names in d (variable name) plus fields label, units, and comment. Dataset m can contain variables not currently being used. To add the labels and units into d and to store comments separately, use the following example.

n <- names(d)
i <- n %nin% m$name
if(any(i)) cat('The following variables have no metadata:',
               paste(n[i], collapse=', '), '\n')
vcomment        <- m$comment
names(vcomment) <- m$name
mn              <- subset(m, name %in% n)
labs            <- mn$label
un              <- mn$units
names(labs)     <- names(un) <- mn$name
d <- upData(d, labels=labs, units=un)

To look up the comment for a variable at any time use e.g. vcomment['height']. All comments were saved in vector vcomment in case the metadata dictionary m defined variables that are not in d but were to be imported later in the script. Comments for multiple variables can be looked up using e.g. vcomment[.q(height, weight, age)].

If you want to look up a variable’s comment without having to quote its name use the following:

vcom <- function(...)
  vcomment[as.character(sys.call()[-1])]
# Example usage: vcom(age,sbp,dbp)

The built-in function in R for reading .csv files is read.csv. The Hmisc package has a function csv.get which calls read.csv but offers some enhancements in variable naming, date handling, and reading variable labels from a specified row number. Illegal characters in variable names are changed to periods, and by default underscores are also changed to periods. If any variable names are changed and the labels argument is not given, original variable names are stored in the variable label attributes.

The fread function in the data.table package is blazing fast for reading large files and offers a number of options. csv.get uses it if data.table is installed on the system.

If reading data exported from REDCap that are placed into the project directory I run the following to get rid of duplicate (factor and non-factor versions of variables REDCap produces) variables and automatically convert dates to Date variables:

require(Hmisc)
getRs('importREDCap.r', put='source')  # source() code to define function
mydata <- importREDCap()  # by default operates on last downloaded export
saveRDS(mydata, 'mydata.rds', compress='xz')

When file names are not given to importREDCap the function looks for the latest created .csv file and .R file with same prefix and uses those. See this for more information.

SAS, Stata, and SPSS binary files are converted to R data.frames using the R haven package. Here’s an example:

require(haven)   # after you've installed the package
d <- read_sas('mydata.sas7bdat')
d <- read_xpt('mydata.xpt')        # SAS transport files
d <- read_dta('mydata.dta')        # Stata
d <- read_sav('mydata.sav')        # SPSS
d <- read_por('mydata.por')        # Older SPSS files

These import functions carry variable labels into the data frame and convert dates and times appropriately. Character vectors are not converted to factors.

One of the most important principles to following in programming data analyses is to not do the same thing more than once. Repetitive code wastes time and is harder to maintain. One example of avoiding repetition is in reading a large number of files in R. If the files are stored in one directory and have a consistent file naming scheme (or you want to import every .csv file in the directory), one can avoid naming the individual files. The results may be stored in an R list that has named elements, and there are many processing tasks that can be automated by looping over this list.

In the following example assume that all the data files are .csv files in the current working directory, and they all have names of the form xz*.csv. Let’s read all of them and put each file into a data frame named by the characters in front of .csv. These data frames are stuffed into a list named X. The Hmisc csv.get function is used to read the files, automatically translating dates to Date variables, and because lowernames=TRUE is specified, variable names are translated to lower case. There is an option to fetch variable labels from a certain row of each .csv file but we are not using that.

files <- list.files(pattern='xz.*.csv')  # vector of qualifying file names
# Get base part of *.csv
dnames <- sub('.csv', '', files)
X <- list()
i <- 0
for(f in files) {
  cat('Reading file', f, '\n')
  i <- i + 1
  d <- csv.get(f, lowernames=TRUE)
  # To read SAS, Stata, SPSS binary files use a haven function instead
  # To convert to data.table do setDT(d) here
  X[[dnames[i]]] <- d
}
saveRDS(X, 'X.rds', compress='xz')  # Efficiently store all datasets together

To process one of the datasets one can do things like summary(X[[3]]) or summary(X$baseline) where the third dataset stored was named baseline because it was imported from baseline.csv.

Now there are many possibilities for processing all the datasets at once such as the following.

k   <- length(X)          # number of datasets
nam <- lapply(X, names)   # list with k elements, each is a vector of names
# Get the union of all variable names used in any dataset
sort(unique(unlist(nam))) # all the variables appearing in any dataset
# Get variable names contained in all datasets
common <- names(X[[1]])
for(i in 2 : k) {
  common <- intersect(common, names(X[[i]]))
  if(! length(common)) break  # intersection already empty
}
sort(common)
# Compute number of variables across datasets
nvar <- sapply(X, length)  # or ncol
# Print number of observations per dataset
sapply(X, nrow)
# For each variable name count the number of datasets containing it
w <- data.table(dsname=rep(names(X), nvar), vname=unlist(nam))
w[, .N, keyby=vname]
# For each variable create a comma-separated list of datasets
# containing it
w[, .(datasets=paste(sort(dsname), collapse=', ')), keyby=vname]
# For datasets having a subject ID variable named id compute
# the number of unique ids
uid <- function(d) if('id' %in% names(d)) length(unique(d$id)) else NA
sapply(X, uid)
# To repeatedly analyze one of the datasets, extract it to a single data frame
d <- X$baseline
describe(d)

The Hmisc package cleanup.import function improves imported data storage in a number of ways including converting double precision variables to integer when originally double but not containing fractional values (this halves the storage requirement). Hmisc::upData is the go-to function for annotating data frames/tables, renaming variables, and dropping variables. Hmisc::dataframeReduced removes problematic variables, e.g., those with a high fraction of missing values or that are binary with very small prevalence.

Variable Naming

I prefer short but descriptive variable names. As exemplified above, I use variable labels and units to provide more information. For example I wouldn’t name the age variable age.at.enrollment.yrs but would name it age with a label of Age at Enrollment and with units of years. Short, clear names unclutter code, especially formulas in statistical models. One can always fetch a variable label while writing a program (e.g., typing label(d$age) at the console) to check that you have the right variable (or put the data dictionary in a window for easy reference, as shown below). Hmisc package graphics and table making functions such as summaryM and summary.formula specially typeset units in a smaller font.

Data Dictionary

The Hmisc package contents function will provide a concise data dictionary. Here is an example using the permanent version of the dataset created above, which can be accessed with the Hmisc getHdata function. The top of the contents output has the number of levels for factor variables hyperlinked. Click on the number to go directly to the list of levels for that variable.

The permanent version of stressEcho coded binary variables as 0/1 instead of N/Y.
require(Hmisc)
getHdata(stressEcho)
d <- stressEcho
html(contents(d), levelType='table')
d Contents

Data frame:d

558 observations and 31 variables, maximum # NAs:0  
NameLabelsUnitsLevelsStorage
bhrBasal heart ratebpminteger
basebpBasal blood pressuremmHginteger
basedpBasal Double Product bhr*basebpbpm*mmHginteger
pkhrPeak heart ratemmHginteger
sbpSystolic blood pressuremmHginteger
dpDouble product pkhr*sbpbpm*mmHginteger
doseDose of dobutamine givenmginteger
maxhrMaximum heart ratebpminteger
pctMphrPercent maximum predicted heart rate achieved%integer
mbpMaximum blood pressuremmHginteger
dpmaxdoDouble product on max dobutamine dosebpm*mmHginteger
dobdoseDobutamine dose at max double productmginteger
ageAgeyearsinteger
gender2integer
baseEFBaseline cardiac ejection fraction%integer
dobEFEjection fraction on dobutamine%integer
chestpainChest paininteger
restwmaResting wall motion abnormality on echocardiograminteger
posSEPositive stress echocardiograminteger
newMINew myocardial infarctioninteger
newPTCARecent angioplastyinteger
newCABGRecent bypass surgeryinteger
deathinteger
hxofHTHistory of hypertensioninteger
hxofDMHistory of diabetesinteger
hxofCigHistory of smoking3integer
hxofMIHistory of myocardial infarctioninteger
hxofPTCAHistory of angioplastyinteger
hxofCABGHistory of coronary artery bypass surgeryinteger
any.eventDeath, newMI, newPTCA, or newCABGinteger
ecgBaseline electrocardiogram diagnosis3integer
VariableLevels
gendermale
female
hxofCigheavy
moderate
non-smoker
ecgnormal
equivocal
MI

You can write the text output of contents into a text file in your current working directory, and click on that file in the RStudio Files window to create a new tab in the editor panel where you can view the data dictionary at any time. This is especially helpful if you need a reminder of variable definitions that are stored in the variable labels. Here is an example where the formatted data dictionary is saved.

Users having the xless system command installed can pop up a contents window at any time by typing xless(contents(d)) in the console. xless is in Hmisc.
capture.output(contents(d), file='contents.txt')

Or put the html version of the data dictionary into a small browser window to which you can refer at any point in analysis coding.

cat(html(contents(d)), file='contents.html')
browseURL('contents.html', browser='vivaldi -new-window')

RStudio provides a nice way to do this, facilitated by the htmlView function which may be fetched from Github using the Hmisc getRs function as shown below. htmlView takes any number of objects for which an html method exists to render them. They are rendered in the RStudio Viewer pane. If you are running outside RStudio, your default browser will be launched instead.

Occasionally RStudio Viewer will drop its arrow button making it impossible to navigate back and forth to different html outputs.
Code for htmlView and htmlViewx may be viewed in reptools.r.
getRs('reptools.r', put='source')
# reptools.r defines htmlView, htmlViewx, kabl, maketabs, dataChk
htmlView(contents(d))

In some cases it is best to have a browser window open to the full descriptive statistics for a data table/frame (see below; the describe function also shows labels, units, and levels).

For either approach it would be easy to have multiple tabs open, one tab for each of a series of data tables, or use htmlView.

To be able to have multiple windows open to see information about datasets it is advantageous to open an external browser window. The htmlViewx function defined when you read htmlView.r will by default open a Vivaldi browser window with the first output put in a new window and all subsequent objects displayed as tabs within that same window. This behavior can be controlled with the tab argument, and you can use a different browser by issuing for example options(vbrowser='firefox'). As an example suppose that two datasets were read from the hbiostat.org/data data repository, and the data dictionary and descriptive statistics for both datasets were to be converted to html and placed in an external browser window for the duration of the R session.

In Windows you may have to specify a full path and firefox.exe. The tab names will not be correct until Hmisc 4.7-1 appears.
getHdata(support)
getHdata(support2)
htmlViewx(contents(support ), describe(support ),
          contents(support2), describe(support2))

A screenshot of the result is here.

Data Checking

Besides useful descriptive statistics exemplified below, it is important to flag suspicious values in an automated way. Since checking multiple columns may involve a large number of R expressions to run to classify observations as suspicious, let’s automate the process somewhat by specifying a vector of expressions. Then we have R “compute on the language” to parse the expressions for finding observations to flag, and for printing. This is done by the dataChk function in Github.

The following code results in separate output for each individual data check, in separate Quarto tabs. The dataset does not have a subject ID variable so let’s create one, and also add a site variable to print. Arguments are specified to dataChk so that no tab is produced for a condition that never occurred in the data, and a tab is produced showing all data flags, sorted by id and site.

getRs('reptools.r', put='source')  # Define dataChk, maketabs
require(data.table)
w <- d
setDT(w)
w[, id := 1 : .N]
set.seed(1)
w[, site := sample(LETTERS[1:6], .N, replace=TRUE)]
checks <- expression(
  age < 30 | age > 90,
  gender == 'female' & maxhr > 170,
  baseEF %between% c(72, 77),
  baseEF > 77,
  baseEF > 99,
  sbp > 250 & maxhr < 160)
dataChk(w, checks, id=c('id', 'site'),
        omit0=TRUE, byid=TRUE, html=TRUE)
     id site age 
 1:  14    B  29 
 2:  23    F  91 
 3:  30    D  26 
 4:  60    A  91 
 5:  64    D  92 
 6: 116    A  28 
 7: 235    F  91 
 8: 259    D  29 
 9: 313    C  93 
 
     id site gender maxhr 
 1:  11    A female   171 
 2:  89    B female   182 
 3: 412    E female   200 
 
     id site baseEF 
 1:  56    B     72 
 2: 200    A     75 
 3: 272    A     72 
 4: 366    E     74 
 5: 406    B     77 
 6: 433    A     74 
 7: 495    D     72 
 8: 496    E     74 
 
     id site baseEF 
 1: 299    D     79 
 2: 434    E     83 
 
     id site sbp maxhr 
 1:  51    B 309   146 
 2: 146    E 283   135 
 3: 353    D 274   117 
 
      id site                            Check     Values 
  1:  11    A gender == "female" & maxhr > 170 female 171 
  2:  14    B              age < 30 | age > 90         29 
  3:  23    F              age < 30 | age > 90         91 
  4:  30    D              age < 30 | age > 90         26 
  5:  51    B          sbp > 250 & maxhr < 160    309 146 
  6:  56    B                  baseEF [72, 77]         72 
  7:  60    A              age < 30 | age > 90         91 
  8:  64    D              age < 30 | age > 90         92 
  9:  89    B gender == "female" & maxhr > 170 female 182 
 10: 116    A              age < 30 | age > 90         28 
 11: 146    E          sbp > 250 & maxhr < 160    283 135 
 12: 200    A                  baseEF [72, 77]         75 
 13: 235    F              age < 30 | age > 90         91 
 14: 259    D              age < 30 | age > 90         29 
 15: 272    A                  baseEF [72, 77]         72 
 16: 299    D                      baseEF > 77         79 
 17: 313    C              age < 30 | age > 90         93 
 18: 353    D          sbp > 250 & maxhr < 160    274 117 
 19: 366    E                  baseEF [72, 77]         74 
 20: 406    B                  baseEF [72, 77]         77 
 21: 412    E gender == "female" & maxhr > 170 female 200 
 22: 433    A                  baseEF [72, 77]         74 
 23: 434    E                      baseEF > 77         83 
 24: 495    D                  baseEF [72, 77]         72 
 25: 496    E                  baseEF [72, 77]         74 
      id site                            Check     Values 
 
                              Check n 
 1              age < 30 | age > 90 9 
 2 gender == "female" & maxhr > 170 3 
 3                  baseEF [72, 77] 8 
 4                      baseEF > 77 2 
 5                      baseEF > 99 0 
 6          sbp > 250 & maxhr < 160 3

Descriptive Statistics

The Hmisc describe function is my main tool for getting initial descriptive statistics and quality controlling the data in a univariate fashion. Here is an example. The Info index is a measure of the information content in a numeric variable relative to the information in a continuous numeric variable with no ties. A very low value of Info will occur when a highly imbalanced variable is binary. Clicking on Glossary on the right will pop up a browser window with a more in-depth glossary of terms used in Hmisc package output. It links to hbiostat.org/R/glossary.html which you can link from your reports that use Hmisc.

Info comes from the approximate formula for the variance of a log odds ratio for a proportional odds model/Wilcoxon test, due to Whitehead.
Glossary

Gmd in the output stands for Gini’s mean difference—the mean absolute difference over all possible pairs of different observations. It is a very interpretable measure of dispersion that is more robust than the standard deviation.

w <- describe(d)
html(w)
d Descriptives
d

31 Variables   558 Observations
bhr: Basal heart rate bpm
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580690.99975.2916.57 54.0 58.0 64.0 74.0 84.0 95.3102.0
lowest : 42 44 45 46 47 , highest: 108 115 116 127 210
basebp: Basal blood pressure mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580940.998135.323.35104.0110.0120.0133.0150.0162.3170.1
lowest : 85 88 90 97 98 , highest: 192 194 195 201 203
basedp: Basal Double Product bhr*basebp bpm*mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
55804411101812813 6607 7200 8400 9792116631361014770
lowest : 5000 5220 5280 5400 5460 , highest: 17604 17710 17748 21082 27300
pkhr: Peak heart rate mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
55801051120.625.36 81.85 90.70106.25122.00135.00147.00155.15
lowest : 52 61 62 63 66 , highest: 170 171 176 182 210
sbp: Systolic blood pressure mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
55801421146.940.72 96102120141170200210
lowest : 40 60 70 79 80 , highest: 240 250 274 283 309
dp: Double product pkhr*sbp bpm*mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580508117634576510256113411403317060206442453626637
lowest : 5100 5940 7490 8100 8360 , highest: 32518 33400 33840 38205 45114
dose: Dose of dobutamine given mg
image
nmissingdistinctInfoMeanGmd
558070.8433.758.334
lowest : 10 15 20 25 30 , highest: 20 25 30 35 40 
 Value         10    15    20    25    30    35    40
 Frequency      2    28    47    56    64    61   300
 Proportion 0.004 0.050 0.084 0.100 0.115 0.109 0.538
maxhr: Maximum heart rate bpm
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
55801031119.424.64 82.0 91.0104.2120.0133.0146.0154.1
lowest : 58 62 63 66 67 , highest: 170 171 176 182 200
pctMphr: Percent maximum predicted heart rate achieved %
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580780.99978.5716.86 53 60 69 78 88 97104
lowest : 38 39 40 41 42 , highest: 116 117 126 132 133
mbp: Maximum blood pressure mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
55801320.99915635.03110.0120.0133.2150.0175.8200.0211.1
lowest : 84 90 92 93 96 , highest: 240 250 274 283 309
dpmaxdo: Double product on max dobutamine dose bpm*mmHg
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580484118550538511346128651526018118212392489327477
lowest : 7130 8100 8360 9240 9280 , highest: 32518 33400 33840 38205 45114
dobdose: Dobutamine dose at max double product mg
image
nmissingdistinctInfoMeanGmd
558080.94130.2410.55
lowest : 5 10 15 20 25 , highest: 20 25 30 35 40 
 Value          5    10    15    20    25    30    35    40
 Frequency      7     7    55    73    71    78    62   205
 Proportion 0.013 0.013 0.099 0.131 0.127 0.140 0.111 0.367
age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580620.99967.3413.4146.8551.0060.0069.0075.0082.0085.00
lowest : 26 28 29 30 33 , highest: 89 90 91 92 93
gender
nmissingdistinct
55802
 Value        male female
 Frequency     220    338
 Proportion  0.394  0.606
baseEF: Baseline cardiac ejection fraction %
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580540.99455.610.7132405257626566
lowest : 20 21 22 23 25 , highest: 74 75 77 79 83
dobEF: Ejection fraction on dobutamine %
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580600.99265.2412.3840.049.762.067.073.076.080.0
lowest : 23 25 26 27 28 , highest: 86 87 89 90 94
chestpain: Chest pain
nmissingdistinctInfoSumMeanGmd
558020.641720.30820.4272
restwma: Resting wall motion abnormality on echocardiogram
nmissingdistinctInfoSumMeanGmd
558020.7452570.46060.4978
posSE: Positive stress echocardiogram
nmissingdistinctInfoSumMeanGmd
558020.5531360.24370.3693
newMI: New myocardial infarction
nmissingdistinctInfoSumMeanGmd
558020.143280.050180.09549
newPTCA: Recent angioplasty
nmissingdistinctInfoSumMeanGmd
558020.138270.048390.09226
newCABG: Recent bypass surgery
nmissingdistinctInfoSumMeanGmd
558020.167330.059140.1115
death
nmissingdistinctInfoSumMeanGmd
558020.123240.043010.08247
hxofHT: History of hypertension
nmissingdistinctInfoSumMeanGmd
558020.6253930.70430.4173
hxofDM: History of diabetes
nmissingdistinctInfoSumMeanGmd
558020.6992060.36920.4666
hxofCig: History of smoking
image
nmissingdistinct
55803
 Value           heavy   moderate non-smoker
 Frequency         122        138        298
 Proportion      0.219      0.247      0.534
hxofMI: History of myocardial infarction
nmissingdistinctInfoSumMeanGmd
558020.5991540.2760.4004
hxofPTCA: History of angioplasty
nmissingdistinctInfoSumMeanGmd
558020.204410.073480.1364
hxofCABG: History of coronary artery bypass surgery
nmissingdistinctInfoSumMeanGmd
558020.399880.15770.2661
any.event: Death, newMI, newPTCA, or newCABG
nmissingdistinctInfoSumMeanGmd
558020.402890.15950.2686
ecg: Baseline electrocardiogram diagnosis
image
nmissingdistinct
55803
 Value         normal equivocal        MI
 Frequency        311       176        71
 Proportion     0.557     0.315     0.127
 
# To create a separate browser window:
cat(html(w), file='desc.html')
browseURL('desc.html', browser='firefox -new-window')

Better, whether using RStudio or not:

htmlView(w, contents(d))  # or htmlView(describe(d1), describe(d2), ...)
# Use htmlViewx to use an external browser window (see above)

There is also a plot method for describe output. It produces two graphics objects: one for categorical variables and one for continuous variables. The default is to use ggplot2 to produce static graphics.

p <- plot(w, bvspace=2.5)
p$Categorical

p$Continuous

To make creation of tabs in Quarto easier, the maketabs function loops through the elements of a named list and outputs each element into a separate Quarto tab. A wide argument is used to expand the width of the output outside the usual margins. An initblank argument creates a first tab that is empty. This allows one to show nothing until one of the other tabs is clicked.

This function was defined by the earlier getRs('reptools.r', put='source') command and its source code may be viewed here.
maketabs(p)

By specifying grType option you can instead get plotly graphics that use hover text to show more information, especially when hovering over the leftmost dot or tick mark for a variable.

options(grType='plotly')
p <- plot(w, bvspace=2.5)
maketabs(p, wide=TRUE)

See this for other Hmisc functions for descriptive graphics and tables, especially for stratified descriptive statistics for categorical variables. The summaryM function prints a tabular summary of a mix of continuous and categorical variables. Here is an example where stratification is by history of myocardial infarction (MI).

The initblank option to the maketabs function creates the initially opened tab, which has no content, thereby leaving it to the user to click one of the other tabs before output is shown.
require(data.table)
setDT(d)   # turn d into a data table
# tables() with no arguments will give a concise summary of all active data tables
w <- d
w[, hxofMI := factor(hxofMI, 0 : 1, c('No history of MI', 'History of MI'))]
vars <- setdiff(names(d), 'hxofMI')
form <- as.formula(paste(paste(vars, collapse='+'), '~ hxofMI'))
print(form)

bhr + basebp + basedp + pkhr + sbp + dp + dose + maxhr + pctMphr + mbp + dpmaxdo + dobdose + age + gender + baseEF + dobEF + chestpain + restwma + posSE + newMI + newPTCA + newCABG + death + hxofHT + hxofDM + hxofCig + hxofPTCA + hxofCABG + any.event + ecg ~ hxofMI

s <- summaryM(form, data=d, test=TRUE)
w <- list('Table 1'                   = html(s, exclude1=TRUE, npct='both', digits=3, middle.bold=TRUE),
          'Categorical Variable Plot' = plot(s, which='categorical', vars=1 : 4, height=600, width=1000),
          'Continuous Variable Plot'  = plot(s, which='continuous',  vars=1 : 4))
# initblank=TRUE: Put an empty tab first so that nothing initially displays
maketabs(w, wide=TRUE, initblank=TRUE)
Descriptive Statistics (N=558).
No history of MI
N=404		 History of MI
N=154		 Test Statistic
Basal heart rate
bpm
 65 74 85 63 72 84 F1 556=1.41, P=0.2351
Basal blood pressure
mmHg
 120 134 150 120 130 150 F1 556=1.39, P=0.2381
Basal Double Product bhr*basebp
bpm*mmHg
 8514 9874 11766 8026 9548 11297 F1 556=3.32, P=0.0691
Peak heart rate
mmHg
 107 123 136 104 120 132 F1 556=2.35, P=0.1261
Systolic blood pressure
mmHg
 124 146 174 115 134 158 F1 556=12.1, P<0.0011
Double product pkhr*sbp
bpm*mmHg
 14520 17783 21116 13198 15539 18885 F1 556=15, P<0.0011
Dose of dobutamine given
mg
: 10		 0.00 2⁄404 0.00 0⁄154 χ26=8.77, P=0.1872
  15 0.05 21⁄404 0.05 7⁄154
  20 0.10 40⁄404 0.05 7⁄154
  25 0.11 45⁄404 0.07 11⁄154
  30 0.11 43⁄404 0.14 21⁄154
  35 0.11 45⁄404 0.10 16⁄154
  40 0.51 208⁄404 0.60 92⁄154
Maximum heart rate
bpm
 107 122 134 102 118 130 F1 556=4.05, P=0.0451
Percent maximum predicted heart rate achieved
%
 69.0 78.0 89.0 70.0 77.0 87.5 F1 556=0.5, P=0.4791
Maximum blood pressure
mmHg
 138 154 180 130 142 162 F1 556=13, P<0.0011
Double product on max dobutamine dose
bpm*mmHg
 15654 18666 21664 14489 16785 19680 F1 556=16.1, P<0.0011
Dobutamine dose at max double product
mg
: 5		 0.01 4⁄404 0.02 3⁄154 χ27=8.5, P=0.292
  10 0.01 6⁄404 0.01 1⁄154
  15 0.11 43⁄404 0.08 12⁄154
  20 0.14 58⁄404 0.10 15⁄154
  25 0.13 54⁄404 0.11 17⁄154
  30 0.13 51⁄404 0.18 27⁄154
  35 0.10 40⁄404 0.14 22⁄154
  40 0.37 148⁄404 0.37 57⁄154
Age
years
 59.0 68.0 75.0 63.2 71.0 76.8 F1 556=9.75, P=0.0021
gender : female 0.68 273⁄404 0.42 65⁄154 χ21=30, P<0.0012
Baseline cardiac ejection fraction
%
 55 59 63 40 54 60 F1 556=56.4, P<0.0011
Ejection fraction on dobutamine
%
 65.0 70.0 74.2 50.0 64.5 70.0 F1 556=50.3, P<0.0011
Chest pain 0.29 119⁄404 0.34 53⁄154 χ21=1.29, P=0.2572
Resting wall motion abnormality on echocardiogram 0.57 230⁄404 0.18 27⁄154 χ21=69.7, P<0.0012
Positive stress echocardiogram 0.21 86⁄404 0.32 50⁄154 χ21=7.56, P=0.0062
New myocardial infarction 0.03 14⁄404 0.09 14⁄154 χ21=7.4, P=0.0072
Recent angioplasty 0.02 10⁄404 0.11 17⁄154 χ21=17.8, P<0.0012
Recent bypass surgery 0.05 21⁄404 0.08 12⁄154 χ21=1.35, P=0.2462
death 0.04 15⁄404 0.06 9⁄154 χ21=1.23, P=0.2672
History of hypertension 0.69 280⁄404 0.73 113⁄154 χ21=0.89, P=0.3462
History of diabetes 0.36 147⁄404 0.38 59⁄154 χ21=0.18, P=0.6742
History of smoking : heavy 0.21 83⁄404 0.25 39⁄154 χ22=3.16, P=0.2062
  moderate 0.24 96⁄404 0.27 42⁄154
  non-smoker 0.56 225⁄404 0.47 73⁄154
History of angioplasty 0.04 15⁄404 0.17 26⁄154 χ21=28.4, P<0.0012
History of coronary artery bypass surgery 0.08 34⁄404 0.35 54⁄154 χ21=59.6, P<0.0012
Death, newMI, newPTCA, or newCABG 0.12 48⁄404 0.27 41⁄154 χ21=18.1, P<0.0012
Baseline electrocardiogram diagnosis : normal 0.59 240⁄404 0.46 71⁄154 χ22=8, P=0.0182
  equivocal 0.29 117⁄404 0.38 59⁄154
  MI 0.12 47⁄404 0.16 24⁄154
a b c represent the lower quartile a, the median b, and the upper quartile c for continuous variables.
Tests used: 1Wilcoxon test; 2Pearson test .

Semi-interactive stratified spike histograms are also useful descriptive plots. These plots also contain a superset of the quantiles used in box plots, and the legend is clickable, allowing any of the statistical summaries to be turned off.

d[, histboxp(x=maxhr, group=ecg, bins=200)]

Data Manipulation and Aggregation

The data.table package provides a concise, consistent syntax for managing simple and complex data manipulation tasks, and it is extremely efficient for large datasets. One of the best organized tutorials is this, and a cheatsheet for data transformation is here. A master cheatsheet for data.table is here from which the general syntax below is taken, where DT represents a data table.

   DT[ i,  j,  by ] # + extra arguments
        |   |   |
        |   |    -------> grouped by what?
        |    -------> what to do?
         ---> on which rows?

Data tables are also data frames so work on any method handling data frames. But data tables do not contain rownames.

Several data.table examples follow. I like to hold the current dataset in d to save typing. Some basic operations on data tables are:

d[2]                  # print 2nd row
d[2:4]                # print rows 2,3,4
d[y > 2 & z > 3]      # rows satisfying conditions
d[, age]              # retrieve one variable
d[, .(age, gender)]   # make a new table with two variables
i <- 2; d[, ..i]      # get column 2
v <- 'age'; d[, ..v]  # get variable named by contents of v
d[.N]                 # last row
d[, .N, keyby=age]    # number of rows for each age, sorted
d[1:10, bhr:pkhr]     # first 10 rows, variables bhr - pkhr
d[1:10, !(bhr:pkhr)]  # all but those variables
d[, 2:4]              # get columns 2-4

Analyzing Selected Variables and Subsets

d[, html(describe(age))]
Descriptives age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580620.99967.3413.4146.8551.0060.0069.0075.0082.0085.00
lowest : 26 28 29 30 33 , highest: 89 90 91 92 93 
d[, html(describe(~ age + gender))]
age + gender Descriptives
age + gender

2 Variables   558 Observations
age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580620.99967.3413.4146.8551.0060.0069.0075.0082.0085.00
lowest : 26 28 29 30 33 , highest: 89 90 91 92 93
gender
nmissingdistinct
55802
 Value        male female
 Frequency     220    338
 Proportion  0.394  0.606
 
d[gender == 'female', html(describe(age))]   # analyze age for females
Descriptives age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
3380580.99967.0113.7447.0050.7059.2568.0076.0083.0085.00
lowest : 26 28 29 33 34 , highest: 87 88 89 90 91 
html(describe(d[, .(age, gender)], 'Age and Gender Stats'))
Age and Gender Stats Descriptives
Age and Gender Stats

2 Variables   558 Observations
age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
5580620.99967.3413.4146.8551.0060.0069.0075.0082.0085.00
lowest : 26 28 29 30 33 , highest: 89 90 91 92 93
gender
nmissingdistinct
55802
 Value        male female
 Frequency     220    338
 Proportion  0.394  0.606
 
# Separate analysis by female, male.  Use keyby instead of by to sort the usual way.
d[, print(describe(age, descript=gender)), by=gender]
male : Age [years] 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     220        0       51    0.999    67.86    12.91    45.95    51.00 
     .25      .50      .75      .90      .95 
   62.00    69.00    75.00    81.00    86.00 

lowest : 30 34 38 40 43, highest: 88 89 91 92 93
female : Age [years] 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     338        0       58    0.999    67.01    13.74    47.00    50.70 
     .25      .50      .75      .90      .95 
   59.25    68.00    76.00    83.00    85.00 

lowest : 26 28 29 33 34, highest: 87 88 89 90 91
Empty data.table (0 rows and 1 cols): gender
# Compute mean and median age by gender
d[, .(Mean=mean(age), Median=median(age)), by=gender]
   gender     Mean Median
1:   male 67.85909     69
2: female 67.00888     68
# To create a new subset
w <- d[gender == 'female' & age < 70, ]

Adding or Changing Variables

With data.table you can create a new data table with added variables, or you can add or redefine variables in a data table in place. The latter has major speed and memory efficiency implications when processing massive data tables. Here d refers to a different data table from the one used above.

# Rename a variable
setnames(d, c('gender', 'height'),
            c(  'sex',      'ht'))
# Easier:
setnames(d, .q(gender, height),
            .q(   sex,     ht))

# Or
ren <- .q(gender=sex, height=ht)
setnames(d, names(ren), ren)

# Or
rename <- function(x, n) setnames(x, names(n), n)
rename(d, .q(gender=sex, height=ht))

# For all variables having baseline_ at the start of their names, remove it
n <- names(d)
setnames(d, n, sub('^baseline_', '', n))   # ^ = at start; use $ for at end

# For all variables having baseline_ at the start of their names, remove it
# and add Baseline to start of variable label
n    <- names(d)
n    <- n[grepl('^baseline_', n)]
ren  <- sub('^baseline_', '', n); names(ren) <- n
# Fetch vector of labels for variables whose names start with baseline_
labs <- sapply(d, label)[n]   # label() result is "" if no label
labs <- paste('Baseline', labs)
d    <- updata(d, rename=ren, labels=labs)

# Change all names to lower case
n <- names(d)
setnames(d, n, tolower(n))

# Abbreviate names of all variables longer than 10 characters
# abbreviate() will ensure that all names are unique
n <- names(d)
setnames(d, n, abbreviate(n, minlength=10))

# For any variables having [...] or (...) in their labels, assume these
# are units of measurement and move them from the label to the units
# attribute
d <- upData(d, moveUnits=TRUE)

# Add two new variables, storing
result in a new data table
z <- d[, .(bmi=wt / ht ^ 2, bsa=0.016667 * sqrt(wt * ht))]

# Add one new variable in place
d[, bmi := wt / ht ^ 2]

# Add two new variables in place
d[, `:=`(bmi = wt / ht ^ 2, bas=0.016667 * sqrt(wt * ht))]
d[, .q(bmi, bsa) := .(wt / ht ^ 2, 0.016667 * sqrt(wt * ht))]

# Compute something requiring a different formula for different types
# of observations
d[, htAdj := ifelse(sex == 'male', ht, ht * 1.07)]  # better" use fifelse in data.table
d[, htAdj := ht * ifelse(sex == 'male', 1, 1.07)]
d[, htAdj := (sex == 'male') * ht + (sex == 'female') * ht * 1.07]
d[sex == 'male',   htAdj := ht]
d[sex == 'female', htAdj := ht * 1.07]
d[, htAdj := fcase(sex == 'male',   ht,          # fcase is in dta.table
                   sex == 'female', ht * 1.07)]
d[, htAdj := fcase(sex = 'female', ht * 1.07, default = ht)]

# Add label & optional units (better to use upData which works on data tables)
adlab <- function(x, lab, un='') {
  label(x) <- lab
  if(un != '') units(x) <- un
  x
}
d[, maxhr := adlab(maxhr, 'Maximum Heart Rate', '/m')]

# Delete a variable (or use upData)
d[, bsa := NULL]

# Delete two variables
d[, `:=`(bsa=NULL, bmi=NULL)]
d[, .q(bsa, bmi) := NULL]

Recoding Variables

# Group levels a1, a2 as a & b1,b2,b3 as b
d[, x := factor(x, .q(a1,a2,b1,b2,b3),
                   .q( a, a, b, b, b))]
# Regroup an existing factor variable
levels(d$x) <- list(a=.q(a1,a2), b=.q(b1,b2,b3))
# or
d <- upData(d, levels=list(x=list(a=.q(a1,a2), b=.q(b1,b2,b3))))
# Or manipulate character strings
d[, x := substring(x, 1, 1)]   # replace x with first character of levels
# or
levels(d$x) <- substring(levels(d$x), 1, 1)

# Recode a numeric variable with values 0, 1, 2, 3, 4 to 0, 1, 1, 1, 2
d[, x := 1 * (x %in% 1:3) + 2 * (x == 4)]
d[, x := fcase(x %in% 1:3, 1,          # fcase is in data.table
               x == 4,     2)]
d[, x := fcase(x %between% c(1,3), 1,
               x    ==       4,    2)]  # %between% is in data.table

# Recode a series of conditions to a factor variable whose value is taken
# from the last condition that is TRUE using Hmisc::score.binary
# Result is a factor variable unless you add retfactor=FALSE
d[, x := score.binary(x1 == 'symptomatic',
                      x2 %in% .q(stroke, MI),
                      death)]
# Same but code with numeric points
d[, x := score.binary(x1 == 'symptomatic',
                      x2 %in% .q(stroke, MI),
                      death,  # TRUE/FALSE or 1/0 variable
                      points=c(1,2,10))]
# Or just reverse the conditions and use fcase which stops at the first
# condition met
d[, x := fcase(death,                  'death',        # takes precedence
               x2 %in% .q(stroke, MI), 'stroke/MI',    # takes next precedence
               x1 == 'symptomatic',    'symptomatic',
               default =               'none')]

# Recode from one set of character strings to another using named vector
# A named vector can be subscripted with character strings as well as integers
states <- c(AL='Alabama', AK='Alaska', AZ='Arizona', ...)
# Could also do:
#  states <- .q(AL=Alabama, AK=Alaska, AZ=Arizona, ..., NM='New Mexico', ...) 
# or do a merge for table lookup (see later)
d[, State := states[state]]
# states are unique, state can have duplicates and all are recoded
d[, State := fcase(state == 'AL', 'Alabama',  state='AK', 'Alaska', ...)]

# Recode from integers 1, 2, ..., to character strings
labs <- .q(elephant, giraffe, dog, cat)
d[, x := labs[x]]

# Recode from character strings to integers
d[, x := match(x, labs)]
d[, x := fcase(x == 'elephant', 1, 
               x == 'giraffe',  2,
               x == 'dog',      3,
               x == 'cat',      4)]

As an example of more complex hierarchical recoding let’s define codes in a nested list.

a <- list(plant =
            list(vegetable = .q(spinach, lettuce, potato),
                 fruit     = .q(apple, orange, banana, pear)), 
          animal =
            list(domestic  = .q(dog, cat, horse),
                 wild      = .q(giraffe, elephant, lion, tiger)) )
a
$plant
$plant$vegetable
[1] "spinach" "lettuce" "potato" 

$plant$fruit
[1] "apple"  "orange" "banana" "pear"  


$animal
$animal$domestic
[1] "dog"   "cat"   "horse"

$animal$wild
[1] "giraffe"  "elephant" "lion"     "tiger"   
a <- unlist(a)
a
plant.vegetable1 plant.vegetable2 plant.vegetable3     plant.fruit1 
       "spinach"        "lettuce"         "potato"          "apple" 
    plant.fruit2     plant.fruit3     plant.fruit4 animal.domestic1 
        "orange"         "banana"           "pear"            "dog" 
animal.domestic2 animal.domestic3     animal.wild1     animal.wild2 
           "cat"          "horse"        "giraffe"       "elephant" 
    animal.wild3     animal.wild4 
          "lion"          "tiger" 
# Pick names of unlist'ed elements apart to define kingdom and type
n <- sub('[0-9]*$', '', names(a))  # remove sequence numbers from ends of names
# Names are of the form kingdom.type; split at .
s       <- strsplit(n, '.', fixed=TRUE) 
kingdom <- sapply(s, function(x) x[1])
type    <- sapply(s, function(x) x[2])
# or:   (note \\. is escaped . meaning not to use as wild card)
#        .* = wild card: any number of characters
# kingdom <- sub('\\..*', '', n)  # in xxx.yyy remove .yyy
# type    <- sub('.*\\.', '', n)  # in xxx.yyy remove xxx.
names(kingdom) <- names(type) <- a
w <- data.table(kingdom, type, item=a, key=c('kingdom', 'item'))
w
    kingdom      type     item
 1:  animal  domestic      cat
 2:  animal  domestic      dog
 3:  animal      wild elephant
 4:  animal      wild  giraffe
 5:  animal  domestic    horse
 6:  animal      wild     lion
 7:  animal      wild    tiger
 8:   plant     fruit    apple
 9:   plant     fruit   banana
10:   plant vegetable  lettuce
11:   plant     fruit   orange
12:   plant     fruit     pear
13:   plant vegetable   potato
14:   plant vegetable  spinach
# Example table lookups
cat(kingdom['dog'], ':', type['dog'], '\n')
animal : domestic 
kingdom[.q(dog, cat, spinach)]
     dog      cat  spinach 
"animal" "animal"  "plant" 
type   [.q(dog, cat, giraffe, spinach)]
        dog         cat     giraffe     spinach 
 "domestic"  "domestic"      "wild" "vegetable" 
# But what if there is a plant named the same as an animal?
# Then look up on two keys
w[.('animal', 'lion'), type]
[1] "wild"

Computing Summary Statistics

Many applications can use the automatically created data.table object .SD which stands for the data table for the current group being processed. If .SDcols were not specified, all variables would be attempted to be analyzed. Specify a vector of variable names as .SDcols to restrict the analysis. If there were no by variable(s), .SD stands for the whole data table.

# Compute the number of distinct values for all variables
nd <- function(x) length(unique(x))
d[, sapply(.SD, nd)]
      bhr    basebp    basedp      pkhr       sbp        dp      dose     maxhr 
       69        94       441       105       142       508         7       103 
  pctMphr       mbp   dpmaxdo   dobdose       age    gender    baseEF     dobEF 
       78       132       484         8        62         2        54        60 
chestpain   restwma     posSE     newMI   newPTCA   newCABG     death    hxofHT 
        2         2         2         2         2         2         2         2 
   hxofDM   hxofCig    hxofMI  hxofPTCA  hxofCABG any.event       ecg 
        2         3         2         2         2         2         3 
# Same but only for variables whose names contain hx and either D or M
d[, sapply(.SD, nd), .SDcols=patterns('hx', 'D|M')]
hxofDM hxofMI 
     2      2 
# Compute means on all numeric variables
mn <- function(x) mean(x, na.rm=TRUE)
d[, lapply(.SD, mn), .SDcols=is.numeric]
        bhr   basebp   basedp     pkhr      sbp       dp     dose    maxhr
1: 75.29032 135.3244 10181.31 120.5502 146.9158 17633.84 33.75448 119.3692
    pctMphr mbp  dpmaxdo  dobdose      age   baseEF    dobEF chestpain
1: 78.56989 156 18549.88 30.24194 67.34409 55.60394 65.24194 0.3082437
     restwma     posSE      newMI   newPTCA    newCABG      death    hxofHT
1: 0.4605735 0.2437276 0.05017921 0.0483871 0.05913978 0.04301075 0.7043011
      hxofDM  hxofPTCA  hxofCABG any.event
1: 0.3691756 0.0734767 0.1577061 0.1594982
# Compute means on all numeric non-binary variables
nnb <- function(x) is.numeric(x) && length(unique(x)) > 2
d[, lapply(.SD, mn), .SDcols=nnb]
        bhr   basebp   basedp     pkhr      sbp       dp     dose    maxhr
1: 75.29032 135.3244 10181.31 120.5502 146.9158 17633.84 33.75448 119.3692
    pctMphr mbp  dpmaxdo  dobdose      age   baseEF    dobEF
1: 78.56989 156 18549.88 30.24194 67.34409 55.60394 65.24194
# Print frequency tables of all categorical variables with > 2 levels
cmult <- function(x) ! is.numeric(x) && length(unique(x)) > 2
tab <- function(x) {
  z <- table(x, useNA='ifany')
  paste(paste0(names(z), ': ', z), collapse=', ')
}
d[, lapply(.SD, tab), .SDcols=cmult]
                                      hxofCig
1: heavy: 122, moderate: 138, non-smoker: 298
                                   ecg
1: normal: 311, equivocal: 176, MI: 71

Tabulate all variables having between 3 and 10 distinct values and create a side effect when data.table is running that makes the summarization function tab store all values and frequencies in a growing list Z so that kable can render a markdown table after we pad columns to the maximum length of all columns (maximum number of distinct values).

# Using <<- makes data.table have a side effect of augmenting Z and
# Align in the global environment
tab <- function(x) {
  z <- table(x, useNA='ifany')
  i <- length(Z)
  Z[[i+1]] <<- names(z)
  Z[[i+2]] <<- as.vector(z)
  Align <<- c(Align, if(is.numeric(x)) 'r' else 'l', 'r')
  length(z)
}
discr <- function(x) { i <- length(unique(x)); i > 2 & i < 11 }
#                                       or i %between% c(2,11)
Z    <- list(); Align <- character(0)
w    <- d[, lapply(.SD, tab), .SDcols=discr]
maxl <- max(w)
# Pad shorter vectors with blanks
Z <- lapply(Z, function(x) c(x, rep('', maxl - length(x))))
Z <- do.call('cbind', Z)  # combine all into columns of a matrix
colnames(Z) <- rep(names(w), each=2)
colnames(Z)[seq(2, ncol(Z), by=2)] <- 'Freq'
knitr::kable(Z, align=Align)
dose Freq dobdose Freq hxofCig Freq ecg Freq
10 2 5 7 heavy 122 normal 311
15 28 10 7 moderate 138 equivocal 176
20 47 15 55 non-smoker 298 MI 71
25 56 20 73
30 64 25 71
35 61 30 78
40 300 35 62
40 205

A better approach is to let the kables function put together a series of separate markdown tables of different sizes. By using the “updating Z in the global environment” side effect we are able to let data.table output any type of objects of non-conformable dimensions over variables (such as frequency tabulations).

tab <- function(x) {
  z <- table(x, useNA='ifany')
  i <- length(Z)
  w <- matrix(cbind(names(z), z), ncol=2,
              dimnames=list(NULL, c(vnames[i+1], 'Freq')))
  Z[[i+1]] <<- knitr::kable(w, align=c(if(is.numeric(x)) 'r' else 'l', 'r'))
  length(z)
}
discr <- function(x) { i <- length(unique(x)); i > 2 & i < 11 }
Z      <- list()
vnames <- names(d[, .SD, .SDcols=discr])
w      <- d[, lapply(.SD, tab), .SDcols=discr]
knitr::kables(Z)
dose Freq
10 2
15 28
20 47
25 56
30 64
35 61
40 300
dobdose Freq
5 7
10 7
15 55
20 73
25 71
30 78
35 62
40 205
hxofCig Freq
heavy 122
moderate 138
non-smoker 298
ecg Freq
normal 311
equivocal 176
MI 71

Use a similar side-effect approach to get separate html describe output by gender.

g <- function(x, by) {
  Z[[length(Z) + 1]] <<- describe(x, descript=paste('age for', by))
  by
}
Z <- list()
by <- d[, g(age, gender), by=gender]
# Make Z look like describe() output for multiple variables
class(Z) <- 'describe'
attr(Z, 'dimensions') <- c(nrow(d), nrow(by))
attr(Z, 'descript') <- 'Age by Gender'
html(Z)
Age by Gender Descriptives
Age by Gender

2 Variables   558 Observations
age for male: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
2200510.99967.8612.9145.9551.0062.0069.0075.0081.0086.00
lowest : 30 34 38 40 43 , highest: 88 89 91 92 93
age for female: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
3380580.99967.0113.7447.0050.7059.2568.0076.0083.0085.00
lowest : 26 28 29 33 34 , highest: 87 88 89 90 91 
# Compute a 1-valued statistic on multiple variables, by cross-classification
# of two variables.  Do this on a subset.  .SDcols=a:b uses variables in order
# Use keyby instead of by to order output the usual way
d[age < 70, lapply(.SD, mean), keyby=.(gender, newMI), .SDcols=pkhr:dp]
   gender newMI     pkhr      sbp       dp
1:   male     0 122.0561 147.7664 17836.24
2:   male     1 115.6667 140.6667 16437.67
3: female     0 122.8492 150.5084 18509.03
4: female     1 123.5714 171.5714 21506.71
# Compute multiple statistics on one variable
# Note: .N is a special variable: count of obs for current group
d[, .(Max=max(bhr), Min=min(bhr), Mean=mean(bhr), N=.N), by=.(gender, newMI)]
   gender newMI Max Min     Mean   N
1:   male     0 127  42 74.15385 208
2:   male     1  89  59 71.25000  12
3: female     0 210  45 75.96894 322
4: female     1 108  65 79.43750  16
# Same result another way
g <- function(x) list(Max=max(x), Min=min(x), Mean=mean(x), N=length(x))
d[, g(bhr), by=.(gender, newMI)]  # if g returned a vector instead, use as.list(g(bhr))
   gender newMI Max Min     Mean   N
1:   male     0 127  42 74.15385 208
2:   male     1  89  59 71.25000  12
3: female     0 210  45 75.96894 322
4: female     1 108  65 79.43750  16
d[, as.list(quantile(bhr)), by=gender]
   gender 0% 25% 50% 75% 100%
1:   male 42  63  72  83  127
2: female 45  65  75  85  210
# Compute mean bhr by quintiles of age using Hmisc::cut2
# Bad statistical practice; use scatterplot + smoother instead
d[, .(Mean=mean(bhr)), keyby=.(fifth=cut2(age, g=5))]
     fifth     Mean
1: [26,60) 78.54545
2: [60,67) 76.96907
3: [67,72) 74.26316
4: [72,78) 75.92793
5: [78,93] 70.03846
# Compute multiple statistics on multiple variables
d[, lapply(.SD, quantile), by=gender, .SDcols=.q(bhr, pkhr, sbp)]
    gender bhr   pkhr    sbp
 1:   male  42  52.00  80.00
 2:   male  63 106.75 120.00
 3:   male  72 123.00 140.00
 4:   male  83 136.00 166.00
 5:   male 127 176.00 250.00
 6: female  45  61.00  40.00
 7: female  65 106.25 122.25
 8: female  75 122.00 141.50
 9: female  85 134.00 170.00
10: female 210 210.00 309.00
# Similar but put percentile number in front of statistic value
# Do only quartiles
g <- function(x) {
  z <- quantile(x, (1:3)/4, na.rm=TRUE)
  paste(format(names(z)), format(round(z)))
}
d[, lapply(.SD, g), by=gender, .SDcols=.q(bhr, pkhr, sbp)]
   gender    bhr    pkhr     sbp
1:   male 25% 63 25% 107 25% 120
2:   male 50% 72 50% 123 50% 140
3:   male 75% 83 75% 136 75% 166
4: female 25% 65 25% 106 25% 122
5: female 50% 75 50% 122 50% 142
6: female 75% 85 75% 134 75% 170
# To have more control over labeling and to have one row per sex:
g <- function(x) {
  s <- sapply(x, quantile, na.rm=TRUE)  # compute quantiles for all variables -> matrix
  h <- as.list(s)  # vectorizes first
  # Cross row names (percentile names) with column (variable) names
  # paste(b, a) puts variable name in front of percentile
  names(h) <- outer(rownames(s), colnames(s), function(a, b) paste(b, a))
  h
}
d[, g(.SD), by=gender, .SDcols=.q(bhr, pkhr, sbp)]
   gender bhr 0% bhr 25% bhr 50% bhr 75% bhr 100% pkhr 0% pkhr 25% pkhr 50%
1:   male     42      63      72      83      127      52   106.75      123
2: female     45      65      75      85      210      61   106.25      122
   pkhr 75% pkhr 100% sbp 0% sbp 25% sbp 50% sbp 75% sbp 100%
1:      136       176     80  120.00   140.0     166      250
2:      134       210     40  122.25   141.5     170      309
# Restrict to variables bhr - basedp in order columns created in data table
d[, g(.SD), by=gender, .SDcols=bhr : basedp]
   gender bhr 0% bhr 25% bhr 50% bhr 75% bhr 100% basebp 0% basebp 25%
1:   male     42      63      72      83      127        85     120.00
2: female     45      65      75      85      210        90     122.25
   basebp 50% basebp 75% basebp 100% basedp 0% basedp 25% basedp 50% basedp 75%
1:        130        145         203      5220    7976.25       9483   11183.50
2:        136        150         201      5000    8562.00      10063   11891.25
   basedp 100%
1:       16704
2:       27300
# Can put ! in front of a sequence of variables to do the opposite

# To add duplicated means to raw data use e.g.
# d[, Mean := mean(x), by=sex]

Summary Statistics With Marginal Summaries

The cube function in the data.table package will compute all possible marginal statistics. When there is only one by variable, the overall statistic is computed in addition to compute stratified values. When a dimension is being marginalized over, the value of the by variable for that dimension will be NA.

mn  <- function(x) as.double(mean(x, na.rm=TRUE))
# as.double ensures consistency of numeric type across groups
Nna <- function(x) sum(! is.na(x))
cube(d, .(Meanbhr = mn(bhr), N = Nna(bhr)), by='gender', id=TRUE)
   grouping gender  Meanbhr   N
1:        0   male 73.99545 220
2:        0 female 76.13314 338
3:        1   <NA> 75.29032 558
cube(d, .(Meanbhr = mn(bhr), N = Nna(bhr)), by=.q(gender, hxofMI), id=TRUE)
   grouping gender           hxofMI  Meanbhr   N
1:        0   male    History of MI 73.22472  89
2:        0 female No history of MI 76.30403 273
3:        0   male No history of MI 74.51908 131
4:        0 female    History of MI 75.41538  65
5:        1   male             <NA> 73.99545 220
6:        1 female             <NA> 76.13314 338
7:        2   <NA>    History of MI 74.14935 154
8:        2   <NA> No history of MI 75.72525 404
9:        3   <NA>             <NA> 75.29032 558
# id=TRUE creates output variable grouping to detail level of marginalization
# It is a binary representation, e.g. if by has 3 variables and a row
# is marginalizing over the first and third variables,
# grouping=binary 101 = 5
# Use groupingsets() to control the marginalizations
# Example: marginalize only one variable at a time
groupingsets(d, .(Meanbhr = mn(bhr), N=Nna(bhr)), by=.q(gender, hxofMI),
             sets=list('gender', 'hxofMI'), id=TRUE)
   grouping gender           hxofMI  Meanbhr   N
1:        1   male             <NA> 73.99545 220
2:        1 female             <NA> 76.13314 338
3:        2   <NA>    History of MI 74.14935 154
4:        2   <NA> No history of MI 75.72525 404

Merging Data

Consider a baseline dataset b and a longitudinal dataset L, with subject ID of id.

For more information see this, this, this, this and this. To merge any number of datasets at once and obtain a printed report of how the merge went, use the Hmisc Merge function.
b <- data.table(id=1:4, age=c(21, 28, 32, 23), key='id')
L <- data.table(id  = c(2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5),
                day = c(1, 2, 3, 1, 2, 3, 4, 1, 2, 1, 2, 3),
                y    =  1 : 12, key='id')
b
   id age
1:  1  21
2:  2  28
3:  3  32
4:  4  23
L
    id day  y
 1:  2   1  1
 2:  2   2  2
 3:  2   3  3
 4:  3   1  4
 5:  3   2  5
 6:  3   3  6
 7:  3   4  7
 8:  4   1  8
 9:  4   2  9
10:  5   1 10
11:  5   2 11
12:  5   3 12
# Merge b and L to look up baseline age and associate it with all follow-ups
b[L, on=.(id)]   # Keeps all ids in L (left inner join)
    id age day  y
 1:  2  28   1  1
 2:  2  28   2  2
 3:  2  28   3  3
 4:  3  32   1  4
 5:  3  32   2  5
 6:  3  32   3  6
 7:  3  32   4  7
 8:  4  23   1  8
 9:  4  23   2  9
10:  5  NA   1 10
11:  5  NA   2 11
12:  5  NA   3 12
L[b, on=.(id)]   # Keeps all ids in b (right inner join)
    id day  y age
 1:  1  NA NA  21
 2:  2   1  1  28
 3:  2   2  2  28
 4:  2   3  3  28
 5:  3   1  4  32
 6:  3   2  5  32
 7:  3   3  6  32
 8:  3   4  7  32
 9:  4   1  8  23
10:  4   2  9  23
L[b, on=.(id), nomatch=NULL]  # Keeps only ids in both b and L (right outer join)
   id day y age
1:  2   1 1  28
2:  2   2 2  28
3:  2   3 3  28
4:  3   1 4  32
5:  3   2 5  32
6:  3   3 6  32
7:  3   4 7  32
8:  4   1 8  23
9:  4   2 9  23
uid <- unique(c(b[, id], L[, id]))
L[b[.(uid), on=.(id)]]         # Keeps ids in either b or c (full outer join)
    id day  y age
 1:  1  NA NA  21
 2:  2   1  1  28
 3:  2   2  2  28
 4:  2   3  3  28
 5:  3   1  4  32
 6:  3   2  5  32
 7:  3   3  6  32
 8:  3   4  7  32
 9:  4   1  8  23
10:  4   2  9  23
11:  5   1 10  NA
12:  5   2 11  NA
13:  5   3 12  NA
merge(b, L, by='id', all=TRUE) # also full outer join; calls merge.data.table
    id age day  y
 1:  1  21  NA NA
 2:  2  28   1  1
 3:  2  28   2  2
 4:  2  28   3  3
 5:  3  32   1  4
 6:  3  32   2  5
 7:  3  32   3  6
 8:  3  32   4  7
 9:  4  23   1  8
10:  4  23   2  9
11:  5  NA   1 10
12:  5  NA   2 11
13:  5  NA   3 12

For very large data tables, giving the data tables keys will speed execution, e.g.:

setkey(d, id)
setkey(d, state, city)

Join/merge can be used for data lookups:

s <- data.table(st=.q(AL, AK, AZ, CA, OK), y=5:1)
stateAbbrevs <- data.table(state=state.abb, State=state.name)
s[stateAbbrevs, , on=.(st=state), nomatch=NULL]
   st y      State
1: AL 5    Alabama
2: AK 4     Alaska
3: AZ 3    Arizona
4: CA 2 California
5: OK 1   Oklahoma

Reshaping Data

To reshape data from long to wide format, take the L data table above as an example.

w <- dcast(L, id ~ day, value.var='y')
w
   id  1  2  3  4
1:  2  1  2  3 NA
2:  3  4  5  6  7
3:  4  8  9 NA NA
4:  5 10 11 12 NA
# Now reverse the procedure
m <- melt(w, id.vars='id', variable.name='day', value.name='y')
m
    id day  y
 1:  2   1  1
 2:  3   1  4
 3:  4   1  8
 4:  5   1 10
 5:  2   2  2
 6:  3   2  5
 7:  4   2  9
 8:  5   2 11
 9:  2   3  3
10:  3   3  6
11:  4   3 NA
12:  5   3 12
13:  2   4 NA
14:  3   4  7
15:  4   4 NA
16:  5   4 NA
setkey(m, id, day)   # sorts
m[! is.na(y)]
    id day  y
 1:  2   1  1
 2:  2   2  2
 3:  2   3  3
 4:  3   1  4
 5:  3   2  5
 6:  3   3  6
 7:  3   4  7
 8:  4   1  8
 9:  4   2  9
10:  5   1 10
11:  5   2 11
12:  5   3 12

Other Manipulations of Longitudinal Data

For the longitudinal data table L carry forward to 4 days the subject’s last observation on y if it was assessed earlier than day 4.

w <- copy(L)   # fresh start with no propagation of changes back to L
# Only needed if will be using := to compute variables in-place and
# you don't want the new variables also added to L
# This is related to data.table doing things by reference instead of 
# making copies.  w <- L does not create new memory for w.
# Compute each day's within-subject record number and last record number
# Feed this directly into a data.table operation to save last records 
# when the last is on a day < 4
u <- w[, .q(seq, maxseq) := .(1 : .N, .N), by=id][seq == maxseq & day < 4,]
# Extra observations to fill out to day 4
u <- u[, .(day = (day + 1) : 4, y = y), by=id]
u
   id day  y
1:  2   4  3
2:  4   3  9
3:  4   4  9
4:  5   4 12
w <- rbind(L, u, fill=TRUE)
setkey(w, id, day)  # sort and index
w
    id day  y
 1:  2   1  1
 2:  2   2  2
 3:  2   3  3
 4:  2   4  3
 5:  3   1  4
 6:  3   2  5
 7:  3   3  6
 8:  3   4  7
 9:  4   1  8
10:  4   2  9
11:  4   3  9
12:  4   4  9
13:  5   1 10
14:  5   2 11
15:  5   3 12
16:  5   4 12

Find the first time at which y >= 3 and at which y >= 7.

day[y >= 3] is read as “the value of day when y >= 3”. It is a standard subscripting operation in R for two parallel vectors day and y. Taking the minimum value of day satisfying the condition gives us the first qualifying day.
L[, .(first3 = min(day[y >= 3]),
      first7 = min(day[y >= 7])), by=id]
   id first3 first7
1:  2      3    Inf
2:  3      1      4
3:  4      1      1
4:  5      1      1

Same but instead of resulting in an infinite value if no observations for a subject meet the condition, make the result NA.

mn <- function(x) if(length(x)) min(x) else as.double(NA)
# as.double needed because day is stored as double precision
# (type contents(L) to see this) and data.table requires
# consistent storage types
L[, .(first3 = mn(day[y >= 3]),
      first7 = mn(day[y >= 7])), by=id]
   id first3 first7
1:  2      3     NA
2:  3      1      4
3:  4      1      1
4:  5      1      1

Add a new variable z and compute the first day at which z is above 0.5 for two days in a row for the subject. Note that the logic below looks for consecutive days for which records exist for a subject. To also require the days to be one day apart add the clause day == shift(day) + 1 after shift(z) > 0.5.

set.seed(1)
w <- copy(L)
w[, z := round(runif(.N), 3)]
u <- copy(w)
u
    id day  y     z
 1:  2   1  1 0.266
 2:  2   2  2 0.372
 3:  2   3  3 0.573
 4:  3   1  4 0.908
 5:  3   2  5 0.202
 6:  3   3  6 0.898
 7:  3   4  7 0.945
 8:  4   1  8 0.661
 9:  4   2  9 0.629
10:  5   1 10 0.062
11:  5   2 11 0.206
12:  5   3 12 0.177
mn <- function(x)
  if(! length(x) || all(is.na(x))) as.double(NA) else min(x, na.rm=TRUE)
u[, consecutive := z > 0.5 & shift(z) > 0.5, by=id][, 
    firstday    := mn(day[consecutive]),     by=id]
u
    id day  y     z consecutive firstday
 1:  2   1  1 0.266       FALSE       NA
 2:  2   2  2 0.372       FALSE       NA
 3:  2   3  3 0.573       FALSE       NA
 4:  3   1  4 0.908          NA        4
 5:  3   2  5 0.202       FALSE        4
 6:  3   3  6 0.898       FALSE        4
 7:  3   4  7 0.945        TRUE        4
 8:  4   1  8 0.661          NA        2
 9:  4   2  9 0.629        TRUE        2
10:  5   1 10 0.062       FALSE       NA
11:  5   2 11 0.206       FALSE       NA
12:  5   3 12 0.177       FALSE       NA

In general, using methods that involve counters makes logic more clear, easier to incrementally debug, and easier to extend the condition to any number of consecutive times. Create a function that computes the number of consecutive TRUE values or ones such that whenever the sequence is interrupted by FALSE or 0 the counting starts over. As before we compute the first day at which two consecutive z values exceed 0.5.

nconsec is modified from code found here.
nconsec <- function(x) x * sequence(rle(x)$lengths)
# rle in base R: run length encoding
# Example:
x <- c(0,0,1,1,0,1,1,0,1,1,1,1)
nconsec(x)
 [1] 0 0 1 2 0 1 2 0 1 2 3 4
# To require that the time gap between measurements must be <= 2 time
# units use the following example
t <- c(1:9, 11, 14, 15)
rbind(t=t, x=x)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
t    1    2    3    4    5    6    7    8    9    11    14    15
x    0    0    1    1    0    1    1    0    1     1     1     1
nconsec(x & t <= shift(t) + 2)
 [1] 0 0 1 2 0 1 2 0 1 2 0 1
u <- copy(w)
# nconsec(z > 0.5) = number of consecutive days (counting current
# day) for which the subject had z > 0.5
u[, firstday := mn(day[nconsec(z > 0.5) == 2]), by=id] 
#               |  |  |                    |
#         minimum  |  |                    |
#                day  |                    |
#             such that                    |
#  it's the 2nd consecutive day with z > 0.5
u
    id day  y     z firstday
 1:  2   1  1 0.266       NA
 2:  2   2  2 0.372       NA
 3:  2   3  3 0.573       NA
 4:  3   1  4 0.908        4
 5:  3   2  5 0.202        4
 6:  3   3  6 0.898        4
 7:  3   4  7 0.945        4
 8:  4   1  8 0.661        2
 9:  4   2  9 0.629        2
10:  5   1 10 0.062       NA
11:  5   2 11 0.206       NA
12:  5   3 12 0.177       NA

Overlap Joins

The foverlaps function in data.table provides an amazingly fast way to do complex overlap joins. Our first example is modified from an example in the help file for foverlaps. An annotation column is added to describe what happened.

d1 <- data.table(w     =.q(a, a, b, b, b),
                 start = c( 5, 10, 1, 25, 50),
                 end   = c(11, 20, 4, 52, 60))
d2 <- data.table(w     =.q(a, a, b),
                 start = c(1, 15,  1),
                 end   = c(4, 18, 55),
                 name  = .q(dog, cat, giraffe),
                 key   = .q(w, start, end))
f <- foverlaps(d1, d2, type="any")
ann <- c('no a overlap with d1 5-11 & d2 interval',
         'a 10-20 overlaps with a 16-18',
         'b 1-4 overlaps with b 1-55',
         'b 25-62 overlaps with b 1-55',
         'b 50-60 overlaps with b 1-55')
f[, annotation := ann]
f
   w start end    name i.start i.end                              annotation
1: a    NA  NA    <NA>       5    11 no a overlap with d1 5-11 & d2 interval
2: a    15  18     cat      10    20           a 10-20 overlaps with a 16-18
3: b     1  55 giraffe       1     4              b 1-4 overlaps with b 1-55
4: b     1  55 giraffe      25    52            b 25-62 overlaps with b 1-55
5: b     1  55 giraffe      50    60            b 50-60 overlaps with b 1-55
# Don't include record for non-match
foverlaps(d1, d2, type='any', nomatch=NULL)
   w start end    name i.start i.end
1: a    15  18     cat      10    20
2: b     1  55 giraffe       1     4
3: b     1  55 giraffe      25    52
4: b     1  55 giraffe      50    60
# Require the d1 interval to be within the d2 interval
foverlaps(d1, d2, type="within")
   w start end    name i.start i.end
1: a    NA  NA    <NA>       5    11
2: a    NA  NA    <NA>      10    20
3: b     1  55 giraffe       1     4
4: b     1  55 giraffe      25    52
5: b    NA  NA    <NA>      50    60
# Require the intervals to have the same starting point
foverlaps(d1, d2, type="start")
   w start end    name i.start i.end
1: a    NA  NA    <NA>       5    11
2: a    NA  NA    <NA>      10    20
3: b     1  55 giraffe       1     4
4: b    NA  NA    <NA>      25    52
5: b    NA  NA    <NA>      50    60

Now consider an example where there is an “events” dataset e with 0 or more rows per subject containing start (s) and end (e) times and a measurement x representing a daily dose of something given to the subject from s to e. The base dataset b has one record per subject with times c and d. Compute the total dose of drug received between c and d for the subject. This is done by finding all records in e for the subject such that the interval [c,d] has any overlap with the interval [s,e]. For each match compute the number of days in the interval [s,e] that are also in [c,d]. This is given by min(e,d) + 1 - max(c,s). Multiply this duration by x to get the total dose given in [c,d]. For multiple records with intervals touching [c,d] add these products.

base   <- data.table(id    = .q(a,b,c), low=10, hi=20)
events <- data.table(id    = .q(a,b,b,b,k),
                     start = c( 8,  7, 12, 19, 99),
                     end   = c( 9,  8, 14, 88, 99),
                     dose  = c(13, 17, 19, 23, 29))
setkey(base,   id, low,   hi)
setkey(events, id, start, end)
w <- foverlaps(base, events,
               by.x = .q(id, low,   hi),
               by.y = .q(id, start, end ),
               type = 'any', mult='all', nomatch=NA)

w[, elapsed := pmin(end, hi) + 1 - pmax(start, low)]
w[, .(total.dose = sum(dose * elapsed, na.rm=TRUE)), by=id]
   id total.dose
1:  a          0
2:  b        103
3:  c          0

Similar things are can be done with non-equi merges. For those you can require exact subject matches but allow inexact matches on other variables. The following example is modified from here. A medication dataset holds the start and end dates for a patient being on a treatment, and a second dataset visit has one record per subject ID per doctor visit. For each visit look up the drug in effect if there was one.

medication <-
  data.table(ID         = c( 1, 1, 2, 3, 3),
             medication = .q(a, b, a, a, b),
             start      = as.Date(c("2003-03-25","2006-04-27","2008-12-05",
                                    "2004-01-03","2005-09-18")),
             stop       = as.Date(c("2006-04-02","2012-02-03","2011-05-03",
                                    "2005-06-30","2010-07-12")),
             key        = 'ID')
medication
   ID medication      start       stop
1:  1          a 2003-03-25 2006-04-02
2:  1          b 2006-04-27 2012-02-03
3:  2          a 2008-12-05 2011-05-03
4:  3          a 2004-01-03 2005-06-30
5:  3          b 2005-09-18 2010-07-12
set.seed(123)
visit <- data.table(
  ID   = rep(1:3, 4),
  date = sample(seq(as.Date('2003-01-01'), as.Date('2013-01-01'), 1), 12),
  sbp  = round(rnorm(12, 120, 15)),
  key  = c('ID', 'date'))
visit
    ID       date sbp
 1:  1 2004-06-09 113
 2:  1 2006-06-06 147
 3:  1 2008-01-16 126
 4:  1 2009-09-28 127
 5:  2 2003-07-14 138
 6:  2 2006-02-15 122
 7:  2 2006-06-21 127
 8:  2 2009-11-15 101
 9:  3 2005-11-03  91
10:  3 2009-02-04 110
11:  3 2011-03-05 125
12:  3 2012-03-24 112
# Variables named in inequalities need to have variables in
# medication listed first
m <- medication[visit, on = .(ID, start <= date, stop > date)]
m
    ID medication      start       stop sbp
 1:  1          a 2004-06-09 2004-06-09 113
 2:  1          b 2006-06-06 2006-06-06 147
 3:  1          b 2008-01-16 2008-01-16 126
 4:  1          b 2009-09-28 2009-09-28 127
 5:  2       <NA> 2003-07-14 2003-07-14 138
 6:  2       <NA> 2006-02-15 2006-02-15 122
 7:  2       <NA> 2006-06-21 2006-06-21 127
 8:  2          a 2009-11-15 2009-11-15 101
 9:  3          b 2005-11-03 2005-11-03  91
10:  3          b 2009-02-04 2009-02-04 110
11:  3       <NA> 2011-03-05 2011-03-05 125
12:  3       <NA> 2012-03-24 2012-03-24 112
# start and stop dates are replaced with actual date of visit
# drop one of them and rename the other
m[, stop := NULL]
setnames(m, 'start', 'date')
m
    ID medication       date sbp
 1:  1          a 2004-06-09 113
 2:  1          b 2006-06-06 147
 3:  1          b 2008-01-16 126
 4:  1          b 2009-09-28 127
 5:  2       <NA> 2003-07-14 138
 6:  2       <NA> 2006-02-15 122
 7:  2       <NA> 2006-06-21 127
 8:  2          a 2009-11-15 101
 9:  3          b 2005-11-03  91
10:  3          b 2009-02-04 110
11:  3       <NA> 2011-03-05 125
12:  3       <NA> 2012-03-24 112

Graphics

I make heavy use of ggplot2, plotly and R base graphics. plotly is used for interactive graphics, and the R plotly package provides an amazing function ggplotly to convert a static ggplot2 graphics object to an interactive plotly one. If the user goes to the trouble of adding labels for graphics entities (usually points, lines, curves, rectangles, and circles) those labels can become hover text in plotly without disturbing anything in static graphics. As shown here you can sense whether an html or pdf report is being produced, and for html all ggplot2 objects can be automatically transformed to plotly.

With ggplotly extra text appears in front of labels, but the result of ggplotly can be run through Hmisc::ggplotlyr to remove this as shown in the example.

Many types of graphs can be created with base graphics, e.g. hist(age, nclass=50) or Ecdf(age) but using ggplot2 for even simple graphics makes it easy to add handle multiple groups on one graph or to create multiple panels for strata using faceting. ggplot2 has excellent default font sizes and axis labeling that works for most sizes of plots.

Here is a prototypical ggplot2 example illustrating many of the features I most often use. Ignore the ggplot2 label attribute if not using plotly. Options are given to the Hmisc label function so that it will retrieve the variable label and units (if present) and format them for axis labels or tables. The formatting takes into account whether html output is being created and plotly is being used.

# Create a vector of formatted labels for all variables in data
# For variables without labels or units use the variable name
# as the label.  If html and plotly are not in effect use R's
# regular plotmath notation to typeset labels/units
nam   <- names(d)
nv    <- length(nam)
vlabs <- structure(character(nv), names=nam)
for(n in nam)
  vlabs[n] <- label(d[[n]], plot=TRUE, html=ishtml, default=n)

# Define substitutes for xlab and ylab that look up our
# constructed labels.
# Could instead directly use xlab(vlabs['age'])
labx <- function(v) xlab(vlabs[[as.character(substitute(v))]])
laby <- function(v) ylab(vlabs[[as.character(substitute(v))]])
g <-
  ggplot(d, aes(x=age, y=bhr, color=gender, label=paste0('dose:', dose))) +
         geom_point() + geom_smooth() +
         guides(color=guide_legend(title='')) +
         theme(legend.position='bottom') +  # not respected by ggplotly
         labs(caption='Scatterplot of age by basal heart rate stratified by sex') +
         labx(age) + laby(bhr)
# or just xlab('Age in years') + ylab('Basal heart rate')
ggplotlyr(g, remove='.*): ')  # removes paste0("dose:", dose): 
# dose is in hover text for each point

For large datasets the Hmisc package has a function ggfreqScatter that makes it easy to see overlapping points by color coding the frequency of points in each small bin. That way scatterplots scale to very large datasets. Here is an example:

html=TRUE was needed because otherwise axis labels are formatted using R’s plotmath and plotly doesn’t like that.
set.seed(1)
x <- round(rnorm(2000), 1)
y <- 2 * (x > 1.5) + round(rnorm(2000), 1)
z <- sample(c('a', 'b'), 2000, replace=TRUE)
label(x) <- 'X Variable'   # could use xlab() &
label(y) <- 'Y Variable'   # ylab() in ggfreqScatter()
g <- ggfreqScatter(x, y, by=z, html=ishtml)
# If variables were inside a data table use
# g <- d[, ggfreqScatter(x, y, by=z, html=ishtml)]
g

Now convert the graphic to plotly if html is in effect otherwise stay with ggplot2 output.

ggp(g)

When you hover the mouse over a point, its frequency pops up.

Many functions in the Hmisc and rms packages produce plotly graphics directly. One of the most unique functions is dotchartpl.

Analysis

Big Picture

For analysis the sky is the limit, but statistical principles should guide every step. Some of the general principles are

  • If there is to be a pivotal analysis there should be a statistical analysis plan (SAP) for this analysis that does not allow for many “statistician degrees of freedom.” The plan should be completed before doing any analysis that might inform analysis choices in a way that would bias the results (e.g., bias the estimate of treatment effect or bias standard errors of effects in a model).
  • All analyses should be completely reproducible. Explicitly state random number seeds if any random processes (bootstrap, simulation, Bayesian posterior sampling) are involved.
  • Exploratory analysis can take place after any needed SAP is completed.
  • Continuous or ordinal variables should never be dichotomized even for purely descriptive exploratory analysis. For example, computing proportions of patients with disease stratified by quintiles of weight will be both inefficient and misleading.
  • If the study is a randomized trial, presenting descriptive statistics stratified by treatment (“Table 1”) is not helpful, and it is more informative to describe the overall distribution of subjects. Even more helpful is to show how all baseline variables relate to the outcome variable.
  • Above the tendency to interchange the roles of independent and dependent variables by presenting a “Table 2” in such a way that stratifies by the outcome. Stratifying (conditioning) on the outcome is placing it in the role of a baseline variable. Instead, show relationships of baseline variables to outcomes as mentioned in the previous point.
  • Nonparametric smoothers are excellent tools for relating individual continuous variables to an outcome.
  • Models are often the best descriptive tools because they can account for multiple variables simultaneously. For example, instead of computing proportions of missing values of a variable Y stratified by age groups and sex, use a binary logistic regression model to relate smooth nonlinear age and sex to the probability Y is missing.

Replacement for Table 1

Analyses should shed light on the unknown and not dwell on the known. In a randomized trial, the distributions of baseline variables are expected to be the same across treatments, and will be the same once \(N\) is large. What is unknown is how the subject characteristics (and treatment) relate to the outcomes under study. While displaying this trend with a nonparametric smoother, one can simultaneously display the marginal distribution of the characteristic using an extended box plot, spike histogram, or rug plot. Here is an example using the Hmisc summaryRc function. Extended box plots show the mean, median, and quantiles that cover 0.25, 0.5, 0.75, and 0.9 of the distribution.

For more examples see this
getHdata(support)
summaryRc(hospdead ~ age + crea + meanbp + wblc,
          ylim=c(.1, .6), data=support, bpplot='top')

Other Examples

Here is an example in which the distribution of continuous variables is shown, stratified by disease group.

bpplotM(age + crea + meanbp + wblc ~ dzgroup,
        data=support, cex.strip=0.4, cex.means=0.3, cex.n=0.45)

This is better done with interactive plots so that one can for example hover over a corner of a box plot and see which quantile that corner represents.

s <- summaryM(age + crea + meanbp + wblc ~ dzgroup,
              data=support)
options(grType='plotly')
plot(s)

To understand the relationship between two continuous variables we may desire to estimate the mean (using ordinary least squares) and median (using quantile regression) as a smooth function of a predictor (using a restricted cubic spline fit).

require(rms)
ddist <- datadist(support); options(datadist='ddist')
f  <- ols(crea ~ rcs(age, 5), data=support)
p1 <- Predict(f)
g  <- Rq(crea ~ rcs(age, 5), tau=0.5, data=support)
p2 <- Predict(g)
p  <- rbind('OLS Mean'=p1, 'QR Median'=p2)
ggplot(p, groups='.set.')

Trends estimated from this plot are more accurate than those estimated from simply stratifying on age, and there were minimal assumptions. In general the fewest assumptions are made when using known unbiased estimators (in this case the sample mean and sample median) computed on overlapping moving windows (here, of age) that advance by a small amount (here 0.5 years) each time. Estimates are made from the 10th smallest to the 10th largest age to avoid extrapolation. For each target age value the window for estimating creatinine is the target age ± 5y. A data.table non-equi join handles the fuzzy matching into target age intervals.

s      <- as.data.table(support)[, .(age, crea)]
arange <- s[, round(c(sort(age)[10], sort(age, decreasing=TRUE)[10]))]
aseq   <- seq(arange[1], arange[2], by=0.5)
# Same as
# arange <- round(ddist$limits[c('Low:prediction', 'High:prediction'), 'age'])
a <- data.table(tage=aseq, key='tage')    # target ages for estimation
a[, .q(lo, hi) := .(tage - 5, tage + 5)]  # define all windows
m <- a[s, on=.(lo <= age, hi >= age)]     # non-equi join
setkey(m, tage)
g <- function(x) list('Sample Mean'   = mean(x, na.rm=TRUE),
                      'Sample Median' = median(x, na.rm=TRUE))
u <- m[, g(crea), by=tage]
w <- melt(u, id.vars='tage',
          measure.variables=c('Sample Mean', 'Sample Median'),
          variable.name='Statistic', value.name='Creatinine')
# Also compute loess estimates
np <- loess(crea ~ age, data=s)
pc <- predict(np, data.frame(age=aseq))
w  <- rbind(w, data.table(tage=aseq, Creatinine=pc, Statistic='loess'))
# Add estimates from OLS and quantile regression
setDT(p)  # make it a data.table
w  <- rbind(w, p[, .(tage=age, Creatinine=yhat, Statistic=.set.)])
ggplot(w, aes(x=tage, y=Creatinine, color=Statistic)) + geom_line() +
  xlab('Age')

Formatting

I take advantage of special formatting for model fit objects from the rms package by using html or latex methods and putting results='asis' in the chunk header to preserve the formatting. Here is an example.

require(rms)
options(prType='html')  # needed to use special formatting (can use prType='latex')
set.seed(1)
n <- 1000
w <- data.table(x=runif(n), x2=rnorm(n), y=sample(0:1, n, replace=TRUE))
dd <- datadist(w); options(datadist='dd')  # rms needs for summaries and plotting
f <- lrm(y ~ x + rcs(x2, 4), data=w)
f

Logistic Regression Model

 lrm(formula = y ~ x + rcs(x2, 4), data = w)
Model Likelihood
Ratio Test 		Discrimination
Indexes 		Rank Discrim.
Indexes
Obs 1000 LR χ2 1.91 R2 0.003 C 0.526
0 492 d.f. 4 R24,1000 0.000 Dxy 0.051
1 508 Pr(>χ2) 0.7516 R24,749.8 0.000 γ 0.051
max |∂log L/∂β| 2×10-14 Brier 0.249 τa 0.026
β S.E. Wald Z Pr(>|Z|)
Intercept   0.1079  0.2976 0.36 0.7169
x   0.1506  0.2197 0.69 0.4931
x2   0.0502  0.2142 0.23 0.8147
x2’  -0.3705  0.6631 -0.56 0.5763
x2’’   1.3624  2.6219 0.52 0.6033 
anova(f)
Wald Statistics for y
χ2 d.f. P
x 0.47 1 0.4931
x2 1.44 3 0.6961
 Nonlinear 0.33 2 0.8499
TOTAL 1.91 4 0.7523 
latex(f)   # show fitted model in math notation
\[{\rm Prob}\{{\rm y}=1\} = \frac{1}{1+\exp(-X\beta)}, {\rm \ \ where} \\ \] \begin{eqnarray*} X\hat{\beta}= & & \\ & & 0.1078981 +0.1505865\:\mathrm{x} \\ & & + 0.0502146 \mathrm{x_{2}}-0.0323237(\mathrm{x_{2}}+1.665752)_{+}^{3}+0.1188625 (\mathrm{x_{2}}+0.4215632)_{+}^{3} \\ & & -0.1070493(\mathrm{x_{2}}-0.3644079)_{+}^{3}+0.02051052(\mathrm{x_{2}}-1.719818)_{+}^{3} \\ \end{eqnarray*} and \((x)_{+}=x\) if \(x > 0\), 0 otherwise

Caching

The workhorse behind Rmarkdown and Quarto is knitr, which processes the code chunks and properly mingles code and tabular and graphical output. knitr has a built-in caching mechanism to make it so that code is not needlessly executed when the code inputs have not changed. This easy-to-use process does have two disadvantages: the dependencies are not transparent, and the stored cache files may be quite large. I like to take control of caching. To that end, the runifChanged function was written. Here is an example of its use. First a function with no arguments must be composed. This is the (usually slow) function that will be conditionally run if any of a group of listed objects has changed since the last time it was run. This function when needed to be run produces an object that is stored in binary form in a user-specified file (the default file name is the name of the current R code chunk with .rds appended).

# Read the source code for the hashCheck and runifChanged functions from
# https://github.com/harrelfe/rscripts/blob/master/hashCheck.r
getRs('hashCheck.r', put='source')
g <- function() {
  # Fit a logistic regression model and bootstrap it 500 times, saving
  # the matrix of bootstrapped coefficients
  f <- lrm(y ~ x1 + x2, x=TRUE, y=TRUE, data=dat)
  bootcov(f, B=500)
}
set.seed(3)
n   <- 2000
dat <- data.table(x1=runif(n), x2=runif(n),
                  y=sample(0:1, n, replace=TRUE))
# runifChanged will write runifch.rds if needed (chunk name.rds)
# Will run if dat or source code for lrm or bootcov change
b <- runifChanged(g, dat, lrm, bootcov)
dim(b$boot.Coef)
[1] 500   3
head(b$boot.Coef)
       Intercept          x1          x2
[1,]  0.02007292 -0.30079958  0.32416398
[2,]  0.06150624 -0.35741054  0.25522669
[3,]  0.25225861 -0.40094541  0.09290729
[4,]  0.13766665 -0.48661991  0.19684403
[5,] -0.22018456  0.02132711  0.33973578
[6,]  0.18217417 -0.36140896 -0.04873320

Parallel Computing

The runParallel function makes it easy to use available multiprocessor cores to speed up parallel computations especially for simulations. By default it runs the number of available cores, less one. runParallel makes the parallel package easier to use and does recombinations over per-core batches. The user writes a function that does the work on one core, and the same function is run on all cores. This function has set arguments and must return a named list. A base random number seed is given, and the seed is set to this plus i for core number i. The total number of repetitions is given, and this most balanced possible number of repetitions is run on each core to sum to the total desired number of iterations. runifChanged is again used, to avoid running the simulations if no inputs have changed.

# Load runParallel function from github
getRs('runParallel.r', put='source')

# Function to do simulations on one core
run1 <- function(reps, showprogress, core) {
  cof <- matrix(NA, nrow=reps, ncol=3,
                dimnames=list(NULL, .q(a, b1, b2)))
  for(i in 1 : reps) {
    y <- sample(0:1, n, replace=TRUE)
    f <- lrm(y ~ X)
    cof[i, ] <- coef(f)
  }
  list(coef=cof)
}
# Debug one core run, with only 3 iterations
n    <- 300
seed <- 3
set.seed(seed)
X    <- cbind(x1=runif(n), x2=runif(n))  # condition on covariates
run1(3)
$coef
              a         b1        b2
[1,] -0.5455330  0.9572896 0.2215974
[2,] -0.2459193  0.3081405 0.1284809
[3,] -0.1391344 -0.2668562 0.1792319
nsim <- 5000
g <- function() runParallel(run1, reps=nsim, seed=seed)
Coefs <- runifChanged(g, X, run1, nsim, seed)
dim(Coefs)
[1] 5000    3
apply(Coefs, 2, mean)
            a            b1            b2 
 0.0020121803 -0.0007277216 -0.0003258610

Simulation

Some of the best ways to validate an analysis are

  • If using any model/feature selection methods use the bootstrap to check whether the selection process is volatile, e.g., your sample size isn’t large enough too support making hard-and-fast selections of predictors/features
  • Use Monte Carlo simulation to check if the correct model or correct predictors are usually selected
  • Simulate a large dataset under a known model and known parameter values and make sure the estimation process you use can recover the true parameter values
  • Simulate the statistical performance of a method under a variety of conditions

Unlike static papers in the literature, simulation can study the performance of a method under conditions that mimic your situation.

When simulating performance of various quantities under various conditions, creating a large number of variables makes the code long and tedious. It is better to to use data frames/tables or arrays to hold everything together. Data frames and arrays also lead to efficient graphics code for summarization.

Data Table Approach

The expand.grid function is useful for generating all combinations of simulation conditions. Suppose we wanted to simulate statistical properties of the maximum absolute value of the sample correlation coefficient from a matrix of all pairwise correlations from truly uncorrelated variables. We do this while varying the sample size n, the number of variables p, and the type of correlation (Pearson’s or Spearman’s, denoted by r and rho). With expand.grid we don’t need a lot of nested for loops. Run 500 simulations for each condition.

nsim <- 500
R <- expand.grid(n=c(10, 20, 50, 100),
                 p=c(2, 5, 10, 20),
                 sim=1 : nsim)
setDT(R)
set.seed(1)
for(i in 1 : nrow(R)) {  # took 4s
  w <- R[i, ]
  n <- w$n
  p <- w$p
  X <- matrix(rnorm(n * p), ncol=p)
  cors    <- cor(X)
  maxr    <- max(abs(cors[row(cors) < col(cors)])) # use upper triangle
  cors    <- cor(X, method='spearman')
  maxrho  <- max(abs(cors[row(cors) < col(cors)]))
  set(R, i, 'maxr',   maxr)    # set is in data.table & is very fast
  set(R, i, 'maxrho', maxrho)  # set will create the variable if needed
  # If not using data.table use this slower approach:
  # R[i, 'maxr'] <- maxr   etc.
}

The simulations could have been cached or parallelized as discussed above.

Compute the mean (over simulations) maximum correlation (over variables) and plot the results.

w <- R[, .(maxr = mean(maxr), maxrho=mean(maxrho)), by=.(n, p)]
# Make data table taller and thinner to put r, rho as different observations
u <- melt(w, id.vars=c('n', 'p'), variable.name='type', value.name='r')
u[, type := substring(type, 4)]   # remove "max"
ps <- c(2, 5, 10, 20)
u[, p := factor(p, ps, paste0('p:', ps))]
g <- ggplot(u, aes(x=n, y=r, col=type)) + geom_jitter(height=0, width=2) +
      ylim(0, 1) +
      facet_wrap(~ p) +
      guides(color=guide_legend(title='')) +
      ylab('Mean Maximum Correlation Coefficient')
plotly::ggplotly(g)

Array Approach

For large problems, storing results in R arrays is more efficient and doesn’t require duplication of values of n and p over simulations. Once the array is created it can be converted into a data table for graphing.

nsim <- 500
ns   <- c(10, 20, 50, 100)
ps   <- c(2, 5, 10, 20)
R <- array(NA, dim=c(nsim, length(ns), length(ps), 2),
               dimnames=list(NULL,
                             n    = as.character(ns),
                             p    = as.character(ps),
                             type = c('r', 'rho')))
dim(R)
[1] 500   4   4   2
dimnames(R)
[[1]]
NULL

$n
[1] "10"  "20"  "50"  "100"

$p
[1] "2"  "5"  "10" "20"

$type
[1] "r"   "rho"
set.seed(1)
for(i in 1 : nsim) {   # took 1s
  for(n in ns) {
    for(p in ps) {
      X <- matrix(rnorm(n * p), ncol=p)
      for(type in c('r', 'rho')) {
        cors <- cor(X, method=switch(type, r = 'pearson', rho = 'spearman'))
        # or           method=fcase(type == 'r', 'pearson', type == 'rho', 'spearman')
        # or           method=fcase(type == 'r', 'pearson', default='spearman')
        # or           method=c(r='pearson', rho='spearman')[type]
        # or           method=.q(r=pearson, rho=spearman)[type]
        # or           method=if(type == 'r') 'pearson' else 'spearman'
        # or           method=ifelse(type == 'r', 'pearson', 'spearman')
        maxr <- max(abs(cors[row(cors) < col(cors)]))
        R[i, as.character(n), as.character(p), type] <- maxr
      }
    }
  }
}
# Compute mean (over simulations) maximum correlation for each condition
m <- apply(R, 2:4, mean)   # preserve dimensions 2,3,4 summarize over 1
# Convert the 3-dimensional array to a tall and thin data table
# Generalizations of row() and col() used for 2-dimensional matrices
# comes in handy: slice.index
dn <- dimnames(m)
u <- data.table(r    = as.vector(m),
                n    = as.numeric(dn[[1]])[as.vector(slice.index(m, 1))],
                p    = as.numeric(dn[[2]])[as.vector(slice.index(m, 2))],
                type = dn[[3]][as.vector(slice.index(m, 3))])
# If doing this a lot you may want to write a dimension expander function
slice <- function(a, i) {
  dn <- all.is.numeric(dimnames(a)[[i]], 'vector')   # all.is.numeric in Hmisc
  dn[as.vector(slice.index(a, i))]
}
u <- data.table(r    = as.vector(m),
                n    = slice(m, 1),
                p    = slice(m, 2),
                type = slice(m, 3))
  
# Plot u using same ggplot code as above

Computing Environment

The following output is created by the command markupSpecs$html$session(), where markupSpecs is defined in the Hmisc package.

 R version 4.2.0 (2022-04-22)
 Platform: x86_64-pc-linux-gnu (64-bit)
 Running under: Pop!_OS 22.04 LTS
 
 Matrix products: default
 BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
 
 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base     
 
 other attached packages:
 [1] rms_6.3-1         SparseM_1.81      plotly_4.10.0     Hmisc_4.7-1      
 [5] ggplot2_3.3.5     Formula_1.2-4     survival_3.3-1    lattice_0.20-45  
 [9] data.table_1.14.2
To cite R in publications use:

R Core Team (2022). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

To cite the Hmisc package in publications use:

Harrell Jr F (2022). Hmisc: Harrell Miscellaneous. R package version 4.7-1, https://hbiostat.org/R/Hmisc/.

To cite the rms package in publications use:

Harrell Jr FE (2022). rms: Regression Modeling Strategies. https://hbiostat.org/R/rms/, https://github.com/harrelfe/rms.

To cite the data.table package in publications use:

Dowle M, Srinivasan A (2021). data.table: Extension of 'data.frame'. R package version 1.14.2, https://CRAN.R-project.org/package=data.table.

To cite the ggplot2 package in publications use:

Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.

To cite the survival package in publications use:

Therneau T (2022). A Package for Survival Analysis in R. R package version 3.3-1, https://CRAN.R-project.org/package=survival.